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Previous  work  on  the  development  and  numerical  implementation  of  the  Bounding 
Surface  Plasticity  model  for  clays  is  discussed.  Modifications  are  made  to  the  hardening 
relationship  to  improve  the  numerical  performance  in  the  tensile  range.  A  rate  equation  for 
the  loading  surface  is  developed.  Modifications  are  made  to  the  invariant  description  of  the 
bounding  surface  to  avoid  numerical  difficulties  in  evaluating  the  derivatives. 

The  Closest  Point  Projection  method  is  described  for  simple  and  general  internal 
variable  plasticity  models.  The  method  is  developed  within  the  classical  plasticity 
framework  and  uses  the  Newton-Raphson  method  to  satisfy  the  implicit  integration  of  the 
rate  equations  and  the  consistency  condition.  An  explicit  treatment  of  the  internal  variables 
is  discussed.  Application  of  this  method  for  the  Bounding  Surface  Plasticity  model  for 
clays  is  developed  by  adding  an  internal  variable  and  using  the  rate  equation  for  the  loading 
surface. 

A  new  algorithm  coined  “the  Reduced  Newton  method”  is  developed  for  the 
Bounding  Surface  Plasticity  model  for  clays.  It  involves  mapping  the  stress  rate  equations, 
internal  variable  rate  equations  and  the  consistency  condition  into  two  nonlinear  equations 
and  integrating  them  with  a  backwards  Euler  formula  using  Newton-Raphson  iteration. 

Comparison  of  predictions  for  a  number  of  sample  problems  is  made  using  the 
trapezoidal  integration.  Closest  Point  and  Reduced  Newton  methods.  “Exact”  solutions  for 
stress  points  that  start  on  the  bounding  surface  are  developed  by  assigning  an  arbitrary 
stress  path  and  calculating  the  corresponding  strains  using  numerical  integration  with  a  tight 
tolerance.  These  “exact”  solutions  are  used  to  evaluate  the  effectiveness  of  the  proposed 
general  numerical  implementation  of  the  Bounding  Surface  model. 

A  standard  effective  stress  interface  is  proposed  for  finite  element  programs  that  use 
the  Newton-Raphson  method.  The  Reduced  Newton  model  is  implemented  within  the 
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DYSAC2  finite  element  program  and  is  used  to  analyze  an  earth  embankment  subjected  to 
earthquake  and  shock  loads. 
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0.  Numerical  Methods  for  Implementing  the  Bounding  Surface 

Plasticity  Model  for  Clays 

In  most  structural  or  geotechnical  engineering  studies,  the  behavior  of  the  object 
being  designed  is  the  focus.  The  object  is  loaded  externally  with  known  and  unknown 
forces  (e.g.,  deadweight,  wind,  soil,  vehicles,  support  reactions,  etc.)  and  the  observed 
behavior  is  in  displacements  (e.g.,  deformation  under  load)  and,  to  some  degree, 
appearance  (e.g.,  concrete  turning  to  rubble).  To  analyze  this  behavior  mathematically,  the 
following  sets  of  equations  must  be  satisfied: 

1)  the  equations  of  equilibrium  (equations  of  motion  for  dynamic  problems), 

2)  the  equations  of  compatibility  (kinematics),  and 

3)  the  constitutive  equations. 

For  quasistatic  problems  the  equilibrium  equations  relate  the  forces  (internal  and 
external)  to  the  stresses.  The  compatibility  equations  provide  the  kinematic  relationship 
between  the  displacements  and  the  strains.  The  constitutive  equations  are  dependent  on  the 
type  of  material  and  relate  the  stresses  to  the  strains.  The  relationships  among  these 
equations  are  shown  in  Figure  0.0-1. 

The  focus  of  this  study  is  on  the  constitutive  equations  and  their  numerical 
implementation  using  the  Finite  Element  Method  (FEM).  It  is  assumed  that  the  finite 
element  program  is  nonlinear,  implicit  and  has  already  addressed  the  equilibrium  and 
compatibility  equations.  The  FEM  is  generally  a  displacement-based  formulation  which 
implies  that  the  constitutive  equations  are  given  strains  and  expected  to  return  stresses  and 
the  tangent  moduli  to  the  global  solution  iteration  procedure. 


The  constitutive  equations  are  developed  to  model  the  behavior  of  the  material  with 
a  given  set  of  measurable  parameters  obtained  from  laboratory  tests.  To  develop  these 
constitutive  equations  the  following  issues  must  be  addressed; 

1)  development  of  the  constitutive  model  to  match  the  material  behavior  under 
consideration,  and 

2)  implementation  of  the  model  within  an  appropriate  solution  algorithm  (usually 
numerical). 

In  the  evolution  of  the  constitutive  model,  a  great  deal  of  effort  is  spent  comparing 
experimental  data  with  numerical  results  to  validate  and  develop  confidence  in  the  model. 
Important  physical  and  geometric  parameters  such  as  void  ratio  for  soils  are  incorporated 
into  the  model.  These  physical  parameters  allow  the  model  to  have  general  application  to  a 
family  of  material  types  (e.g.,  clays  and  silts).  The  constitutive  model  used  in  this  study  is 
the  Bounding  Surface  Plasticity  model  for  clays. 

The  Bounding  Surface  Plasticity  concept  was  developed  at  the  University  of 
California  at  Berkeley  in  the  mid  1970’s  and  it  promised  to  be  computationally  efficient. 
The  concept  provided  for  a  gradual  transition  from  elastic  to  plastic  material  behavior  and 
was  originally  developed  for  metals.  It  became  especially  useful  for  materials  that  have  no 
distinct  yield  points  (e.g.,  soils,  concrete,  etc.).  Bounding  Surface  Plasticity  was  applied 
to  clay  soils  at  the  University  of  California  at  Davis  during  the  mid  1980’s  and  was 
validated  using  traditional  laboratory  and  geotechnical  centrifuge  soils  tests. 

The  subject  for  this  study  is  the  numerical  implementation  of  the  single  ellipse. 
Bounding  Surface  model  for  geotechnical  analysis  of  clay  soils.  The  implementation 


addresses  two  issues: 


1)  efficiency,  and 

2)  accuracy. 

The  importance  of  efficiency  lies  in  the  structure  of  the  solution  algorithm  of  a 
typical  nonlinear,  implicit  finite  element  program  as  illustrated  by  the  nested  loops  shown  in 
Figure  0.0-2.  Evaluation  of  the  constitutive  equations  occurs  within  the  innermost  loop. 
Thus,  efficiency  of  the  constitutive  model  evaluation  significantly  impacts  the  performance 
of  the  global  analysis  of  a  geotechnical  structure. 

Accurate  tangent  moduli  can  enhance  the  performance  of  the  global  solution 
algorithm  by  helping  to  provide  an  optimal  direction  (when  multiplied  by  the  residual 
vector)  for  the  iteration  which  will  improve  the  global  convergence.  Also,  implied  with  the 
issue  of  accuracy  is  robustness.  Even  for  calculations  that  involve  small  solution  time  steps 
in  the  outer  loop,  the  global  iteration  method  can  sometimes  generate  large  trial  strains 
during  the  iteration  process  even  though  the  eventual  incremental  strains  may  be  small. 
These  large  strains  are  provided  as  input  to  the  constitutive  model  and  reasonable  stresses 
and  tangent  moduli  are  expected  to  be  returned  from  the  material  properties  algorithm.  A 
constitutive  model  and  its  numerical  implementation  that  cannot  provide  reasonable  values 
for  large  strain  increments  may  be  useless  even  for  calculations  that  involve  small  time 
steps. 

One  approach  for  dealing  with  the  accuracy/robustness  issue  is  the  use  of  a  uniform 
substepping  method  at  the  material  model  level.  Given  a  set  of  incremental  strains  from  the 
global  iteration,  the  predicted  stresses  and  tangent  moduli  can  be  compared  using  different 
numbers  of  substeps  across  the  increment.  If  the  stresses  from  one  substep  level  are  within 
a  given  tolerance  of  the  previous  level  (e.g.,  one  step  versus  two  substeps)  then  the  finer 


solution  is  accepted.  In  the  global  analysis,  this  assures  that  the  stresses  and  tangent 
moduli  of  neighboring  elements  are  obtained  within  the  same  degree  of  accuracy  in  stress 
space.  The  uniform  substepping  approach  also  has  the  added  advantage  of  seeking  the 
correct  answer  when  Newton-Raphson-based  methods,  which  cannot  distinguish  the 
correct  solution  from  extraneous  solutions  of  a  nonlinear  problem. 

The  original  numerical  implementation  of  the  Bounding  Surface  model  uses 
trapezoidal  integration,  which  is  second  order  accurate  and  relatively  stable.  It  requires 
information  at  the  beginning  and  at  the  end  of  the  step  and  thus  needs  to  iterate.  The 
implementation  also  includes  substepping  to  enhance  accuracy  and  robustness.  Difficulties 
arise  in  certain  analyses  where  the  consistency  condition  is  not  exactly  satisfied  at  the  end 
of  the  increment,  that  is,  where  the  computed  stress  point  falls  outside  the  bounding  surface 
in  stress  space.  This  study  looks  at  alternative  numerical  implementations  that  would 
prevent  this  behavior  and  promote  greater  efficiency  and  accuracy. 

0.1  Report  Layout 

Section  1  briefly  describes  the  theoretical  aspects  of  the  single  ellipse  Bounding 
Surface  model  for  clays.  It  also  describes  theoretical  and  numerical  modifications  that  were 
made  to  improve  the  original  model  and  develops  relationships  that  are  used  in  its  numerical 
implementation. 

Section  2  describes  the  Closest  Point  Projection  algorithm  and  its  application  to  the 
implementation  of  the  Bounding  Surface  model.  The  Closest  Point  method  is  a  general 
three  dimensional  methodology  for  implementing  constitutive  models  and  was  originally 
developed  in  nonlinear  optimization  theory. 

The  Reduced  Newton  algorithm  is  described  in  Section  3.  It  is  a  specific  Newton- 
Raphson-based  method  applied  to  the  Bounding  Surface  model  and  reduces  the  number  of 
differential  equations  to  be  solved  at  the  innermost  iteration  level. 


\  Section  4  compares  both  the  Closest  Point  and  Reduced  Newton  methods  as  well  as 
the  original  trapezoidal  implementation  to  numerically  “exact”  solutions  of  the  constitutive 
equations.  Exact  solutions  include  stress  paths  that  both  start  on  and  within  the  bounding 
surface  and  highlights  the  behavior  of  the  solution  methods. 

Section  5  describes  the  implementation  of  the  Reduced  Newton  model  into  the 
DYSAC2  finite  element  computer  program.  Implementation  issues  are  discussed.  The 
resulting  code  is  then  applied  to  the  solution  of  a  realistic  geotechnical  engineering  problem. 

Conclusions  and  recommendations  are  given  in  Section  6. 


Figure  0.0- 1 .  Relationship  of  Engineering  Mechanics  Equations. 
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Figure  0.0-2.  Nesting  Levels  of  a  Typical  Implicit  Finite  Element  Program. 


1.  Bounding  Surface  Plasticity  Model  for  Clays 

The  Bounding  Surface  plasticity  concept  was  first  introduced  for  modeling  metals  at 
the  University  of  California  at  Berkeley  (UCB)  in  the  mid  1970’s  [Dafalias  and  Popov, 
1975].  The  concept  was  applied  to  clays  at  the  University  of  California  at  Davis  (UCD)  in 
the  early  1980’s  [Dafalias  and  Herrmann,  1980a;  Dafalias  et.  al.  1980b;  Dafalias  et.  al., 
1980c;  Dafalias  1980d;  Dafalias  and  Herrmann,  1980e;  Dafalias  et.  al.,  1981;  Herrmann, 
et.  al.,  1982;  DeNatale,  et.  al.  1983;  Herrmann,  et.  al.  1983a;  Herrmann  and  Mish,  1983b; 
Herrmann,  et.  al.  1983c;  Herrmann  and  Mish,  1983d;  Herrmann  and  Mish,  1983e; 
Herrmann,  et.  al.,  1985;  Kaliakin,  1985;  Shen,  et.  al.,  1986;  Dafalias  and  Herrmann, 

1986;  Herrmann,  et.  al.,  1987].  A  single  ellipse  model  with  an  associated  plastic  flow  was 
developed  for  clays  in  the  mid  1980’s  [Kaliakin,  1985;  Herrmann,  et.  al.;  1985;  Herrmann, 
et.  al.,  1987].  The  model  was  validated  using  traditional  soils  tests  and  centrifuge  tests 
conducted  at  UCD  [Herrmann,  et.  al.,  1982;  Herrmann,  et.  al.,  1987] 

The  following  sections  describe  the  Bounding  Surface  model  and  the  plasticity 
framework  in  which  it  resides.  Section  1.1  briefly  lays  out  the  classical  plasticity 
framework  for  associative  flow.  Section  1.2  discusses  the  Bounding  Surface  Plasticity 
concept.  Section  1.3  describes  the  nonlinear  elastic  volumetric  relationship.  Section  1.4 
presents  the  single  ellipse  Bounding  Surface,  Section  1 .5  details  the  hardening  function  and 
Section  1.6  describes  the  plastic  modulus.  This  work  was  developed  at  UCD  for  clays 
[Kaliakin,  1985]  and  is  provided  for  reference.  A  small  modification  is  made  in  the 
hardening  function  to  avoid  numerical  difficulties  in  Section  1.5.  A  new  relationship  for 
the  loading  surface  is  given  in  Section  1.7  that  lays  out  the  Bounding  Surface  model  in  a 
more  conventional  classical  plasticity  framework.  Finally,  a  modification  to  improve  the 
invariant  form  description  of  the  model  [Kaliakin,  1985]  is  described  in  Section  1.8. 


1.1  Classical  Plasticity  Framework  for  Associative  Plasticity 

Before  introducing  the  Bounding  Surface  Plasticity  model,  the  classical  plasticity 
framework  using  the  associative  flow  rule  [Dafalias,  1990]  is  presented  in  this  section. 
Classical  plasticity  assumes  that  the  material  has  a  region  of  elastic  behavior  within  which 
loading  followed  by  unloading  returns  to  the  original  state.  During  loading,  however,  if 
the  stress  reaches  a  defined  stress  state,  known  as  the  yield  surface,  the  material  begins  to 
yield  or  permanently  deform.  This  yield  surface  in  one  dimension  consists  of  two  points 
and  is  shown  on  a  simple  stress-strain  diagram  in  Figure  1.1-1.  Unloading  after  yielding 
does  not  return  to  the  initial  state,  but  leaves  a  permanent  deformation  and  possibly  residual 
stresses.  The  process  of  yielding  can  also  redefine  the  size,  shape  and/or  location  of  the 
yield  surface,  as  shown  in  Figure  1.1-1.  In  this  case,  the  change  in  the  surface  is  a  result 
of  the  yielding  and  is  a  function  of  the  stress  and/or  strain  history.  The  history  of  the 
change  in  the  surface  at  a  particular  material  point  is  described  by  the  values  of  one  or  more 
internal  variables.  The  yield  surface  can  also  be  described  by  a  number  of  external 
variables  that  are  functions  of  temperature,  humidity,  age  and  other  conditions.  The 
mathematical  framework  for  this  classical  plasticity  model  will  be  described  below. 

The  yield  function  for  a  classical  plasticity  model  is  shown  as  a  one-dimensional 
surface  in  Figure  1.1-1.  For  more  complex  analyses  (two-  and  three-dimensional)  the  yield 
function  is  described  in  stress  space  as  a  multi-dimensional  surface  (see  Figure  1.1-2).  The 
internal  variables  not  only  describe  how  the  surface  grows,  but  how  it  translates,  changes 
shape  and/or  orientation.  The  yield  surface  is  a  function  of  both  the  stresses  and  the 


internal  variables,  and  is  given  as: 


f{cj,q)  =  0  (1.1-1) 

where  /=  yield  function 

a  =  vector  of  stresses,  [cr,; ,  <7^ ,  CT^; ,  ,  •  •  •] 
q  =  vector  of  internal  variables,  ‘ 

Definition  of  the  yield  surface  places  some  restrictions  on  the  stresses.  Stress  states 
can  lie  within  the  surface  (elastic  region)  or  on  the  surface  itself  (plastic  or,  if  unloading, 
elastic).  They  cannot,  however,  exist  outside  of  the  surface.  This  can  be  easily  seen  from 
Figure  1.1-1. 

The  strains  are  assumed  to  decompose  into  elastic  and  plastic  portions.  This 
kinematic  assumption  describing  the  additive  nature  of  the  elastic  and  plastic  strain  rates  is 
given  as: 

£  =  £"+£’’  (1.1-2) 
where  £  =  vector  of  total  strain  rates ,  , . . .] 

=  vector  of  elastic  strain  rates 
£’’  =  vector  of  plastic  strain  rates. 

The  strains  are  described  in  terms  of  rate  equations  (i.e.,  they  evolve  over  time) 
because  time  incorporates  the  history  of  the  plasticity.  Time  can  also  allow  for  a  viscous 
behavior  although  this  is  not  included  in  this  study.  The  stress  rate  is  a  function  of  the 


elastic  strain  rates  only  and  is  given  as: 


<j  =  C  e"  =  c(e-£'’) 

where  &  =  vector  of  stress  rates 


C  =  elastic  constitutive  matrix, 


C. 


(1.1-3) 


Equation  1.1-3  describes  how  the  stress  rates  are  proportional  to  the  elastic  strain 
rates.  As  in  the  one  dimensional  case,  elastic  loading  and  then  unloading  will  return  to  the 
same  state.  When  the  stresses  reach  the  yield  state  (Equation  1.1-1),  plastic  flow 
commences  (i.e.,  plastic  straining  begins).  Plastic  flow  is  described  by  the  increase  in  the 
plastic  strains  and  is  given  by  the  following  rate  equation: 

e'"  =  y  N((j,q)  (1.1-4) 

where  y=  plasticity  parameter;  y>  0,  for  plastic  flow  (f  =  0),  and 


7=  0,  for  elastic  loading  or  unloading  (f<0) 


N((T,  q)  =  2i  function  describing  the  plastic  strain  directions. 


The  associative  flow  rule  assumes  that  the  direction  of  plastic  flow  is  perpendicular 
to  the  yield  surface  which  is  described  in  stress  space  (Equation  1.1-1).  The  rate  equations 
for  the  plastic  strains  are  therefore  a  function  of  the  yield  surface  and  are  given  as: 


£’’  =Y 


I. 

da 


(1.1-5) 


where  —  =  vector  of  yield  function  derivatives 
do 


i.e.. 


da/do’ 
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The  change  in  the  internal  variables  can  be  described  in  terms  of  rate  equations 
which  are  given  as: 

q  =  Yh{a,q)  (1-1-6) 

where  h(G,q)  =  direction  of  internal  variable  rates. 

Because  of  the  restriction  that  the  stress  state  cannot  exist  outside  the  yield  surface 
(Equation  1.1-1),  during  plastic  flow  the  stress  point  must  remain  on  the  surface.  This  is 
described  mathematically  as  the  consistency  condition  and  is  given  as: 

/  =  0  (1.1-7) 

Applying  the  consistency  condition  to  the  yield  state  (Equation  1.1-1)  and  applying 
the  chain  rule  results  in: 


^a  +  ^q  =  0  (1.1-8) 

dG  dq 

Substituting  in  Equation  1.1-6  into  Equation  1.1-8  results  in  and  expression  for  the 
plasticity  parameter  in  terms  of  the  stress  rates  and  is  given  as: 


1 


(1-1-9) 


Combining  Equations  1.1-3,  1.1-5  and  1.1-10  results  in  the  expression  for  the 


plastic  constitutive  matrix.  This  is  given  as: 


<J  = 


C+H{y) 


c 

fill 

ydGj 

'il' 

\ 

c 

/ 

f— 

\  ( 
c 

}  V 

— 1 
.do  J 

£ 


(1.1-11) 


where  H{ y)  =  Heavyside  function;  H( y)  =  0,  7<  0  and  H(y)  =  I,  y>  0. 


In  classical  plasticity,  a  plastic  step  occurs  (for  a  stress  point  starting  on  the  yield 
surface)  when  a  strain  increment  produces  an  elastic  stress  state  outside  of  the  yield 
surface.  This  is  mathematically  given  by  a  positive  plasticity  parameter  (Equation  1.1-10). 
The  plastic  step  requires  integration  of  Equation  1.1-11  (note  that  the  rate  equations  for  the 
plastic  strains  and  internal  variables  were  used  in  the  derivation  of  this  equation). 
Numerical  integration  techniques,  such  as  the  backward  and  forward  Euler  and  the 
trapezoidal  methods,  provide  approximations  for  these  integrations.  The  consistency 
condition  (i.e.,  enforcing  the  stress  point  to  remain  on  the  surface)  is  satisfied  indirectly 
through  the  plasticity  parameter  and  an  accurate  integration.  Since  the  integration  is  an 
approximation,  the  resulting  integration  error  often  fails  to  satisfy  the  consistency  condition 
exactly. 

1.2  Bounding  Surface  Plasticity  Concept 

The  Bounding  Surface  Plasticity  concept  was  introduced  at  the  University  of 
California,  Berkeley  in  the  1970’s  [Dafalias  and  Popov,  1975].  The  motivation  for  the 
concept  was  the  observation  that  for  most  materials  any  stress-strain  curve  (including 
reversals)  eventually  converged  to  well  defined  “bounds”  in  stress-strain  space.  These 
bounds  cannot  be  crossed  but  can  change  position  during  loading.  An  additional 
observation  was  that  the  rate  of  convergence  of  the  stress-strain  curve  to  the  bound 
appeared  to  be  a  function  of  the  distance  of  the  current  state  from  the  bound.  This  concept 


can  be  described  with  a  typical  uniaxial  stress-plastic  strain  response  shown  in  Figure  1.2- 
1.  As  the  stress  approaches  the  bound,  its  rate  of  convergence  or  its  uniaxial  plastic 
modulus  {W)  decreases  until  it  becomes  tangent  with  the  bound.  The  modulus  is  therefore 

a  function  of  the  distance  between  the  current  stress  state  (C7)  and  the  “image”  stress  ( a )  on 
the  bound. 

The  bounding  surface  in  multiaxial  stress  space  is  described  in  a  similar  manner  to  a 
yield  surface  in  “classical”  plasticity  (i.e.,  stress  points  can  exist  on  and  within,  but  not 
outside  the  bounding  surface).  The  unique  feature  of  the  bounding  surface  concept  is  that 
there  is  a  gradual  transition  from  elasticity  to  plasticity  (i.e.,  plasticity  can  occur  within  the 
surface),  unlike  traditional  yield  surface  models  where  plasticity  occurs  only  when  the 
stress  point  is  on  the  surface.  To  accomplish  this,  the  theory  incorporates  particular 
features  that  allow  it  to  operate  within  the  classical  incremental  plasticity  framework.  These 
features  will  be  discussed  in  the  following  paragraphs. 

Bounding  Surface  Plasticity  was  developed  at  UCB  [Dafalias  and  Popov,  1975]  as 
a  means  for  introducing  a  gradual  transition  from  elasticity  to  plasticity.  The  bounding 
surface  is  defined  in  stress  space  and  is  given  as: 

F{a,q)  =  0  (1.2-1) 

where  F  =  bounding  surface  function 
a  =  vector  of  “image”  stresses. 

The  bounding  surface  can  also  have  an  elastic  nucleus  defined  within  it.  The 
current  stress  point  defines  another  surface  known  as  the  loading  surface.  These  surfaces 
are  shown  in  Figure  1.2-2. 

For  “classical”  associative  plasticity,  the  plastic  strain  rates  are  a  function  of  the 
plasticity  parameter  and  are  perpendicular  to  the  yield  surface  at  the  stress  point  (Equation 
1.1-5).  For  stress  points  within  the  bounding  surface  this  information  is  defined  by 


drawing  a  line  from  a  projection  center  (a)  through  the  current  stress  point  and  projecting 

this  to  a  point  on  the  bounding  surface  (see  Figure  1.2-2).  This  is  known  as  the  mapping 
rule.  The  normal  for  the  plastic  strain  rate  is  taken  at  this  “image”  stress  point  on  the 
bounding  surface.  The  strain  rate  is  defined  as; 

?F 

£'’-7^  (1.2-2) 

da 

The  mapping  rule  defines  how  the  current  stress  is  related  to  the  “image”  stress.  It 
is  defined  as: 

a  =  h{a-a)^a  (1-2-3) 

where  b  =  measure  of  distance  between  stress  point  and  surface 


a  =  projection  center. 


Note  that  b  ranges  from  1  (when  the  current  and  “image”  stress  point  coincide)  to  oo 

(at  the  projection  center).  This  parameter  also  implicitly  defines  the  location  of  the  loading 
surface.  The  evolution  of  the  loading  surface  is  discussed  in  Section  1.7. 

The  plasticity  parameter  is  defined  by  applying  the  consistency  condition  to  the 
image  stress  points  on  the  surface.  For  points  within  the  bounding  surface  an  equivalent 
relationship  is  defined  [Kaliakin,  1985]  and  is  given  as: 


1 


(dFA 

— O’ 

Ua  ) 


K 


\  (dF 
- —a 

\da  J 


(1-2-4) 


where 


—  h  =  plastic  modulus  at  the  image  stress  point 
dq 


a  =  image  stress  rate 

Kp  =  plastic  modulus  at  the  current  stress  point 
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<j  =  current  stress  rate. 

Equation  1.2-2  is  the  crux  of  Bounding  Surface  Plasticity.  It  describes  how  the 
plasticity  (via  the  hardening  modulus)  occurs  within  the  bound.  The  plastic  modulus  at  the 

image  stress  point  is  defined  via  the  consistency  condition,  as  in  classical  plasticity, 
with  Equation  1.1-9.  The  plastic  modulus  at  the  current  stress  point  (K^,  within  the 
bounding  surface,  is  a  function  of  (i.e.,  for  stress  points  on  the  bound)  and 

the  distance  to  the  bounding  surface  {b).  This  is  defined  in  more  detail  in  Section  1.6. 

The  following  sections  describe  various  aspects  and  modifications  made  to  the 
Bounding  Surface  Plasticity  concept  for  clays  [Kaliakin,  1985;  Herrmann,  et.  al.,  1985]. 

1.3  Development  of  Nonlinear  Elastic  Volumetric  Relationship 

The  relationship  of  the  elastic  change  in  the  void  ratio  (e)  and  the  volumetric  stress 
(I)  for  unloading-reloading  (URL)  is  modeled  as  a  straight  line  in  log-linear  space 

[Kaliakin,  1985].  This  is  shown  as  the  jcline  in  Figure  1.3-1.  Because  of  problems 

associated  with  I  near  zero  in  log  space,  a  limiting  value  of  the  volumetric  stress  (/, )  is 
given  where  the  relationship  is  changed  from  log  to  linear.  This  is  given  as: 

dl  ^  (/-/;)  +  /, 
de‘  K 

where  I  =  volumetric  stress  (CT;  4-  cr^  +  a^)  (i.e.,  the  first  stress  invariant) 

de^  =  elastic  change  in  the  void  ratio 
I,  =  limiting  value  of  I 

K=  slope  of  unloading  -  reloading  line. 

The  Macaulay  brackets  (  (  )  )  imply  that  (n)  -  n  if  n>0  and  («)  =  0  if  n  <  0. 


This  defines  two  separate  mathematical  regions: 

1)  the  region  of  log-linear  relationship  between  the  void  ratio  and  the  volumetric 
stress  (7  >  /,),  and 

2)  the  region  of  a  linear-linear  relationship  (7  < 

For  this  section,  equations  where  I  >  I,  ait  denoted  as  the  “a”  equations  and  I  <  I,  ait 
denoted  as  the  ”b“  equations.  The  elastic  void  ratio  differentials  for  the  respective  regions 
are  defined  as: 

de^=--dl  (/>/,)  (l-3-2a) 

7 

de‘=--dl  (!<!,)  (1.3-2b) 

h 

The  definition  of  specific  volume  (v)  is  defined  as: 

V  -  l  +  e  (1.3-3) 

Differentiating  Equation  1.3-3  results  in: 

dv  =  de  (1-3-4) 

Equations  1.3-2  now  can  be  written  as  the  elastic  increment  in  specific  volume: 

dv‘=--dl  (l-3-5a) 

7 

dv‘=-  —  dl  (1.3-5b) 

h 


where  dv"=  elastic  change  in  specific  volume. 


The  increment  in  volumetric  strain  can  be  expressed  in  terms  of  the  specific  volume 


for  small  strains  as: 

dd‘=-—  (1-3-6) 

Vo 

where  dd  ‘  =  elastic  change  in  volumetric  strain 
Vg  =  initial  specific  volume. 

Substituting  the  increments  of  specific  volume  (Equations  1.3-5)  into  Equation  1.3- 
6,  the  elastic  increment  of  volumetric  strain  is  given  in  terms  of  the  volumetric  stress 
increment: 

de‘=  —  dl  (1.3-7a) 

(1.3-7b) 

These  relationships  can  be  integrated  from  times  4  to  to  give  the  volumetric 
stress  (7)  in  terms  of  the  elastic  volumetric  strain  {dy. 

(1.3-8a) 

+  (l-3-8b) 

where  p  =  — 

K 

1.4  Single  Ellipse  Bounding  Function  for  Clays 

The  Bounding  Surface  plasticity  concept  was  developed  for  clays  at  the  University 
of  California  at  Davis  (UCD)  in  the  early  1980’s  [Dafalias  and  Herrmann,  1980a;  Dafalias 


et.  al.  1980b;  Dafalias  et.  al.,  1980c;  Dafalias  1980d;  Dafalias  and  Herrmann,  1980e; 
Dafalias  et.  al.,  1981;  Herrmann,  et.  al.,  1982;  DeNatale,  et.  al.  1983;  Herrmann,  et.  al. 
1983a;  Herrmann  and  Mish,  1983b;  Herrmann,  et.  al.  1983c;  Herrmann  and  Mish,  1983d; 
Herrmann  and  Mish,  1983e;  Herrmann,  et.  al.,  1985;  Kaliakin,  1985;  Shen,  et.  al.,  1986; 
Dafalias  and  Herrmann,  1986;  Herrmann,  et.  al.,  1987].  A  single  ellipse  model  with 
associated  plastic  flow  (Section  1.6)  was  developed  for  clays  in  the  mid  1980’s  [Kaliakin, 
1985;  Herrmann,  et.  al.,  1985;  Herrmann,  et.  al.,  1987].  The  model  was  validated  using 
traditional  soils  tests  and  centrifuge  tests  conducted  at  UCD. 

The  single  ellipse  bounding  surface  function  is  expressed  in  terms  of  stress 
invariants  and  is  given  as: 


F=f+{R-iy 


yNj 


R  ”  R 


where  F  =  bounding  function 


(1.4-1) 


I  =b(I-IJ  +  I, 


I  =  first  stress  invariant  (Gj,  +  022  +  ^3s) 


L  =  Ch 

C  =  material  constant  defining  the  projection  center  location 
4  =  bound  size  (i.e.,  intersection  of  bound  with  volumetric  axis) 

J  =bJ 


J  =  second  stress  invariant  = 


V  2 


N  = 


2N. 


l  +  A,,-(l-A,Jsin(3a) 
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a  =  lode  angle  =  jsin 


5  =  third  stress  invariant  = 


7V^  =  slope  of  critical  state  line  in  extension  (7-7  space) 

=  slope  of  critical  state  line  in  compression  (7-7  space) 

7?  =  a  parameter  defining  the  shape  of  the  ellipsoid. 

The  bounding  surface  for  clays  is  shown  in  Figure  1 .4-1 

1.5  Development  of  the  Bounding  Surface  Hardening  Relationship 

The  form  of  the  relationship  for  the  normal  consolidation  line  [Kaliakin,  1985]  is 
similar  to  the  elastic  relationship  (Equation  1.3-1)  and  is  given  as; 

^^JhzlhllL  (1.5-1) 

de  A 

where  7„  =  bounding  surface  size 
e  -  void  ratio 

X  -  slope  of  normal  consolidation  line. 

This  is  shown  in  Figure  1.3-1  as  the  X  line.  The  change  in  the  void  ratio  can  be 
decomposed  into  an  elastic  and  plastic  portion: 

de  =  de'"  +  de^  (1-5-2) 


¥0 


where  d^  -  plastic  change  in  void  ratio. 
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Combining  Equation  1.5-1  with  Equation  1.3-1  (note  that  Equation  1.3-1  is  valid 
for  I  =  I„)  results  in  a  relationship  between  the  bounding  surface  size  and  the  plastic  void 
ratio  and  is: 

^4  _  (A)  ~  (1.5-3) 

de’’  X-K 

In  previous  work  [Kaliakin,  1985]  the  plastic  void  ratio  is  used  as  the  internal 
variable.  The  relationship  describing  its  evolution  is  based  upon  an  associative  flow  rule 
[Kaliakin,  1985]  and  is  given  as: 

&'’=-3v„r^  (1.5-4) 

where  v„  =  initial  specific  volume  =  (1  +  e„) 

7  =  volumetric  image  stress. 

Equation  1.2-1  can  be  rewritten  in  terms  of  volumetric  stress  and  strain  as: 

=  (1-5-5) 

di 

Combining  Equations  1.5-3,  1.5-4  and  1.5-5  results  in  an  expression  relating  the 
plastic  volumetric  strains  to  the  bounding  surface  size.  This  is  given  in  differential  form  as: 

(1-5-6) 

A  —  K 

Once  again,  the  Macaulay  brackets  indicate  two  possible  integrations,  where  I  >  I, 
and  I  <  I,.  However,  when  finding  exact  solutions  (i.e.,  given  a  change  in  and 
calculating  the  plastic  volumetric  strains,  see  Appendix  C),  it  was  noted  that  when  Iq  <  I, 
the  magnitude  of  the  tensile  volumetric  strains  could  not  exceed  a  given  value.  This  is  a 
result  of  the  mathematics  and  not  an  observed  phenomena.  Therefore,  for  this  study  the 
Macaulay  brackets  are  eliminated  in  Equation  1 .5-6  (this  is  consistent  with  the  Closest 
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Point  Projection  application  to  the  Cam-Clay  model  [Simo  and  Meschke,  1993]).  (While 
this  step  may  seem  arbitrary  or  approximate,  it  must  be  remembered  that  the  introduction  of 
was  itself  quite  arbitrary).  The  differential  of  the  bounding  surface  is  now  defined  as: 

(1.5-7) 


where  t  =  — — . 

^  X-K 


Integrating  Equation  1.5-7  between  times  and  results  in  an  expression  for  the 
bound  size  in  terms  of  the  plastic  volumetric  strain: 


(1.5-8) 


1.6  Plastic  Modulus  for  Single  Ellipse  Bounding  Function 


In  classical  plasticity  the  plastic  modulus  is  determined  for  a  stress  point  on  the 
yield  surface  via  the  consistency  condition  (Equation  1.2-3).  This  is  facilitated  by  the  fact 
that  the  stress  point  is  on  the  surface  and  the  yield  function  derivatives  can  be  defined 
directly  by  enforcing  the  consistency  condition.  In  Bounding  Surface  plasticity  the  plastic 
modulus  must  be  defined  within  the  bounding  surface  and  can  not  be  directly  obtained. 

The  behavior  of  the  modulus  can  be  described  by  considering  a  line  from  the  projection 
center  to  the  bounding  surface.  The  modulus  is  a  function  of  the  distance  from  the  bound. 
It  becomes  infinite  at  the  elastic  nucleus  (and  within)  and  approaches  the  classical  plastic 
modulus  definition  at  the  bounding  surface.  The  general  expression  for  the  plastic  modulus 

(A^p)on  the  bound,  at  the  “image”  stress  point,  is  given  by  Equation  1.2-4  and  can  be 
expressed  in  terms  of  the  internal  variable  Ig  as: 


K  =  h 

dl. 


(1.6-1) 
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The  variable  h  is  defined  by  the  rate  equation  of  the  internal  variable,  (see  Equation 
1 . 1-6).  Substitute  Equation  1 .5-5  into  1 .5-7  and  noting  the  form  of  the  rate  equation  for 
the  internal  variable  (Equation  1.1-6),  h  is  given  as: 


h  =  (1.6-2) 

Combining  Equations  1.6-1  and  1.6-2,  the  plastic  modulus  at  the  “image”  stress 
can  be  expressed  as: 


dl  dlf^ 


(1.6-3) 


The  plastic  modulus  {K^  for  points  inside  the  bounding  surface  was  developed  for 
the  single  ellipse  form  of  the  Bounding  Surface  model  for  clays  by  [Kaliakin,  1985].  The 

modulus  is  a  function  of  the  plastic  modulus  at  the  image  stress  point  (isTp)  and  the  distance 

to  the  bound  {b).  It  is  defined  as: 


Jj 


yl  a  + 


sign(nj 


+  1 


dJ 


ip-l) 


(b-{b-i)s;) 


where  =  atmospheric  pressure 


z  = 


JR 

NI. 


m  =  positive  material  constant 

h^{a)  = - 

^  l  +  T]-{l-ri)sm{3a) 


(1.6-4) 


K-K  (f )  ~  value  of  J\  in  triaxial  compression 
=  value  of  in  triaxial  extension 

/ip  =  p-direction  component  of  unit  outward  normal  in  triax  space 
=  shape  hardening  parameter  for  states  near  the  volumetric  axis 

Sp  =  parameter  defining  the  elastic  nucleus  >  l) . 

On  the  bounding  surface  (i.e.,  ^  =  1)  the  plastic  modulus  is  the  same  as  the  image 
modulus  (i.e.,  =  K^).  As  one  moves  toward  the  projection  center  (i.e.,  — >  °°)  the 

value  of  the  term  in  the  Macaulay  brackets  ((  ))  goes  through  zero  and  then  negative  (if 

>  1).  Recalling  the  convention  for  the  Macaulay  brackets,  the  plastic  modulus  becomes 
infinite  as  b  approaches  the  projection  center  or  the  edge  of  elastic  nucleus.  The  elastic 
nucleus  is  defined  by  and  is  given  in  terms  of  the  distance  to  the  bound  by; 


b. 


elastic 


(1.6-5) 


where  =  defines  the  edge  of  the  elastic  nucleus. 

For  stress  points  within  the  elastic  nucleus,  use  of  the  Macaulay  brackets  assures  that  the 
plastic  modulus  remains  at  infinity. 

This  form  of  the  plastic  modulus  was  defined  to  capture  overconsolidated  clay 
behavior  and  to  give  appropriate  behavior  near  the  critical  state  line.  For  further  discussion 
on  the  development  of  this  expression,  the  reader  is  directed  to  previous  work  on  the  single 
ellipse  Bounding  Surface  model  [Kaliakin,  1985]. 


1.7  Development  of  Rate  Equation  for  Loading  Surface 


The  loading  surface  is  the  surface  on  which  the  stress  point  resides  inside  the 
bounding  surface.  The  loading  surface  has  not  been  explicitly  defined  in  the  bounding 


surface  literature  but  is  required  for  any  return  method  algorithm.  The  development  of  the 
rate  equation  for  the  loading  surface  begins  with  the  definition  of  the  image  stress  in 
invariant  space; 
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7  =  Z7(/-/J  +  /^  (1.7-1) 

where 

C  =  projection  center  constant 

b  =  measure  of  distance  from  the  loading  surface  to  the  bounding 
surface. 


J  =  bJ 

a  =  a 

Differentiating  the  image  stress  invariants  with  respect  to  time  yields: 
l^b(i-i,)+b{i-i,)+i, 

J  =  b(J)+b(j) 

a  =  a 


(1.7-2) 

(1.7-3) 

(1.7-4) 

(1.7-5) 

(1.7-6) 


where  7^  = 


The  plasticity  parameter  (Equation  1.2-4)  can  be  expressed  in  terms  of  the 
invariants  as: 


1  (dF  -  dF  ^  dF  ^ 


—  I  +  —  J  +  —OC 
V  dl  dJ  da  J 


I  fdF  ■  dF  ■  I  dF 

- 1  +  ^  7-1--—  a 

\dl  dJ  b  da  J 


(1.7-7) 


Substituting  Equations  1.7-4,  1.7-5  and  1.7-6  into  the  first  part  of  Equation  1.7-7 


yields: 


U/  dJ  )  V 


dF  ■  dF  ■  1  dF  . 

^I  +  r^J  +  -  —  CC 
ol  oJ  b  ooc  j 


dF  ■ 

+  (\-b)^I,  (1.7-8) 


The  terms  in  the  second  parentheses  (multiplied  by  b)  is  recognized  to  be  the 
plasticity  parameter  times  the  hardening  modulus  (i.e.,  the  second  term  of  Equation  1.7-7) 
which  simplifies  Equation  1.7-8  to: 


-  i(dF,  r\  zr  n 


(1.7-9) 


Recall  that  the  projection  center  was  defined  in  Equation  1.7-1  as  =  CIq. 
Combining  Equations  1.5-5  and  1.5-7  (written  in  rate  form)  with  Equation  1.7-1  results  in 
the  expression  for  the  rate  of  the  projection  center  parameter: 


=  3  7  C  Vo 


h  dF 

X-K  dl 


(1.7-10) 


Substituting  Equation  1.7-10  into  1.7-9  the  rate  of  the  loading  surface  {b)  is 
obtained  as: 


K,-bK  -3Cv^  {l-b) 


b  =  Y- 


X-K 


Jfj 


^  dl  dJ 


(1.7-11) 


1.8  Improved  Invariant  Form 

The  model  uses  the  associative  flow  rule  in  which  the  strain  direction  is  given  as  the 
perpendicular  to  the  bounding  surface  (i.e.,  the  first  derivative  of  the  bounding  function 
with  respect  to  the  image  stress).  Since  the  bounding  function  is  given  in  terms  of  stress 


invariants,  this  derivative  is  given  via  the  application  of  the  chain  rule; 


dF  dF  dJ  dF  dJ  dF  da 

— = - ZI - ^  ^  — IT - ^ - IT" 

dl  dcj^j  dJ  dcjy  da  dO^j 


The  lode  angle  (a)  was  used  for  the  third  invariant  in  previous  work  on  Bounding 
Surface  plasticity  [Herrmann,  et.  al.,  1985;  Kaliakin,  1985]  and  Cam-Clay  [ABAQUS, 
1994]  because  of  its  easy  visualization  in  the  7t-plane.  The  partial  derivatives  of  the 
bounding  function  (F)  with  respect  to  the  invariants  can  be  expressed  as: 


31  R  ’ 


(1.8-2) 


1)^  ^ 


(1.8-3) 


da  dN  da  N 


(l- A^^Jcos(3a) 
— 


(1.8-4) 


The  partial  derivatives  of  the  invariants  with  respect  to  the  image  stresses  are  given 


as: 


dl 

do, 


(1.8-5) 


dJ  _ 

2J 


da  V3 


da,  27cos(3a) 


_3^  ^ 

^ik^kj  2  j2 


-H 


7t  ^  ^7t 

—  <a<  — 
6  6 


(1.8-6) 

(1.8-7) 


where 
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Numerical  difficulties  occur  at  the  endpoints  of  the  lode  angle  (i.e.,  a  =  ±f) 


because  the  term,  cos(5a)  =  0.  The  value  of  the  derivatives  associated  with  the  lode  angle 
at  the  endpoints  are: 


—  =  0  (1.8-8) 
da 


da 

da,j 


CO 


(1.8-9) 


The  problem  lies  in  the  fact  that  a  is  discontinuous  at  these  endpoints.  This  can  be 

resolved  for  the  first  derivative  by  reahzing  that  the  cosine  term  cancels  out  when  Equations 
1.8-8  and  1.8-9  are  multiplied  together.  Second  derivative  terms  (as  required  by  the 
Closest  Point  algorithm),  however,  do  not  have  this  obvious  fix.  Another  approach  to  this 

problem  is  noting  that  a  is  not  used  directly  in  the  functional  but  is  used  within 

trigonometric  functions  which  are  continuous.  A  new  variable  is  defined  as: 

jU  =  sin(3a)  (1.8-10) 

Equation  1.8-1  can  now  be  written  in  terms  of  this  new  variable  as: 

dF  dF  dl  dF  dJ  dF  djl 

- =  — —  _  -I — —  _  H - ^ 

dOy  dl  dOy  dJ  dOy  djl  day 


(1.8-11) 


The  derivatives  associated  with  the  third  invariant  (Equations  1.8-4  and  1.8-7)  can 


now  be  written  as: 


72  f 


dll  ^  ’  N 


V  y 


(1.8-12) 


dll  _ -Is 
2/ 


d(7^ 


^ik^kj  2  j2 


45- 

3  ''ij 


(1.8-13) 


This  form  eliminates  the  cos(3a)  term  in  the  first  derivatives  and  thus  avoids 


division  by  zero.  It  also  has  the  added  benefit  of  eliminating  the  cos(3a)  term  in  the  second 


derivatives  (see  Appendix  A). 


First  Effective  Stress  Invariant,  7 


Figure  1.3-1.  Void  Ratio-Volumetric  Stress  Relationship. 


Projection  Center  ' 


Figure  1.4-1.  Single  Ellipse  Bounding  Surface  Model  for  Clays 


2.  Closest  Point  Projection  Method 

The  Closest  Point  Projection  method  was  introduced  to  computational  plasticity  by 
the  late  Professor  Juan  Simo  of  Stanford  University  in  the  early  1990’s  [Simo  and  Hughes, 
1990,  Simo  and  Meschke,  1993].  The  method  was  developed  within  the  classical  plasticity 
framework  and  uses  the  Newton-Raphson  method  in  the  implicit  integration  of  the  rate 
equations  and  the  satisfaction  of  the  consistency  condition.  As  with  any  Newton-Raphson 
method,  a  guess  in  the  “neighborhood”  generally  results  in  quadratic  convergence.  The 
method  has  proven  to  be  extremely  fast  in  12  Flow  models  [Simo  and  Hughes,  1990]. 

In  order  to  employ  the  Newton-Raphson  method,  the  algorithm  requires  that  the 
yield  function  be  twice  differentiable  and  the  internal  variable  rate  equations  be 
differentiable.  Creation  of  this  local  “Jacobian”  requires  additional  calculations,  but 
improves  the  convergence.  The  local  “Jacobian”  also  provides  the  framework  for 
developing  the  algorithmically  consistent  tangent  moduli  [Simo  and  Taylor,  1985]  that 
improves  global  convergence. 

To  remain  consistent  with  Professor  Simo’s  work,  the  differential  form  of  the 

equations  is  used  instead  of  the  rate  form  ( dl  instead  of  /)  throughout  this  section. 

Sections  2.1  and  2.2  lay  out  the  framework  for  the  stress  point  algorithm  and  the  consistent 
moduli,  respectively,  for  a  simple  elastic-plastic  model.  Sections  2.3  and  2.4  describe  the 
stress  point  algorithm  and  the  consistent  moduli  for  a  general  plastic  model  with  internal 
variables.  Section  2.5  extends  the  general  model  using  explicit-implicit  integration. 

Finally,  in  Section  2.6,  application  is  made  of  the  Closest  Point  method  to  the  Bounding 
Surface  Clay  model. 

2.1  Derivation  of  Closest  Point  Stress  Point  Algorithm  (simple  model) 

To  lay  out  the  framework  for  the  Closest  Point  method,  a  simple  elastic-perfectly- 
plastic  model  with  an  associative  flow  rule  will  be  used.  This  will  provide  insight  and 


facilitate  the  application  to  more  complex  models  such  as  the  Bounding  Surface  Clay 
model. 

The  yield  function  for  the  simple  model  is  a  function  of  the  stresses  only  and  is 
given  as: 

f  =  f{a)  =  0  (2.1-1) 

where  /=  yield  function 
(7=  stress  vector. 

The  associative  flow  rule  assumes  that  upon  plastic  flow  the  strains  flow  in  a 
direction  perpendicular  to  the  yield  surface  [Dafalias,  1990]  which  is  described  in  stress 
space.  Thus,  the  rate  equations  for  the  plastic  strains  are  a  function  of  the  yield  surface  (see 
Equation  1.1-5)  and  are  given  as: 

e'’=r^  (2-1-2) 

do 

where  7=  plasticity  parameter 


-2-  =  vector  of  yield  function  derivatives. 
do 

The  desired  result  of  a  constitutive  algorithm  is  the  stresses  at  the  end  of  the  step 
(<^/.+;)  given  a  set  of  strains  (e„+y).  The  stresses  are  a  function  of  only  the  elastic  strains  for 

this  simple  model  and  are  written  as: 


^n+l  —  ^  ^n+l  ^ 


(2.1-3) 


where  (7=  vector  of  stresses  (i.e.,  [  o^,  o^  o,,  ...]^). 


Beside  the  stresses  at  time  t^^,,  the  plastic  strains  are  unknown.  In  order  to 
determine  these,  the  rate  equations  (Equations  2.1-2)  are  integrated  between  times  t„  and 


using  a  backwards  Euler  approximation.  Backwards  Euler  was  chosen  because  it 
involves  the  end  point  values  (i.e.,  at  t^^,)  of  the  various  terms  (e.g.,  plasticity  parameter, 
gradients,  etc.).  This  facilitates  a  simple  derivation  of  the  Closest  Point  Projection  method. 
The  approximation  for  the  plastic  strains  is  given  as: 


e'’  =£^ 

*^n  +  l 


■Jn.X 


K. 

da 


/n+1 


(2.1-4) 


Note  that  will  be  used  for  clarity  throughout  this  section  instead  of  the  more 
appropriate  At 


By  adding  and  subtracting  the  plastic  strains  at  the  beginning  of  the  step  ( )  to  the 


strain  terms  in  the  parenthesis  of  Equation  2.1-3  and  substituting  in  the  right  hand  side  of 
Equation  2.1-4,  the  stress  equations  are  given  in  terms  of  a  “trial”  stress  and  the  unknown 
plasticity  parameter.  The  stress  equations  are  given  as: 


\daj 


(2.1-5) 


where  a“' =  c(£„, -f.'). 

Since  the  total  strains  at  time  are  given  and  the  plastic  strains  at  time  4  are  known  from 
the  previous  calculation,  the  “trial”  stresses  are  in  terms  of  known  quantities  and  can  be 
viewed  as  an  elastic  prediction.  The  concept  of  the  Closest  Point  Projection  method  is  to 
return  the  stresses  back  to  the  surface  via  the  “closest  projection.”  This  is  shown  in  Figure 
2.1-1. 

Since  the  backwards  Euler  method  needs  terms  evaluated  at  time  it  requires 
iteration.  The  Closest  Point  Projection  method  incorporates  Newton-Raphson  iteration  not 
only  to  solve  the  stress  equations  (Equation  2.1-3),  but  to  enforce  the  consistency  condition 


directly.  This  is  accomplished  by  first  defining  local  residuals  using  the  stress  equations. 
The  local  residuals  are  defined  as; 
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/■  ^ 


\daj 


=  0 


n  +  1 


(2.1-6) 


These  residuals  are  differentiated  at  time  (implying  that  values  at  time  r„  are 
constants)  to  give  an  infinitesimal  increment  of  the  residual.  Note  that  the  strains  are 
known  at  times  and  and  therefore  do  not  vary.  Recall,  also,  that  the  yield  equation  in 
this  simple  model  (Equation  2.1-1)  is  a  function  of  stresses  only,  and  the  plasticity 
parameter  is  also  a  variable.  The  resulting  differentiation  is  given  as: 


=dc>„,,+dr,,,C 


£ 


+  7n^^C 


n+] 


do 


(2.1-7) 


The  yield  function  (Equation  2. 1- 1)  can  also  be  differentiated  at  time  By 
setting  this  equation  to  zero,  a  discrete  form  of  the  consistency  condition  (  /  =  0)  is 
established  [Simo  and  Hughes,  1990].  This  is  equivalent  to  satisfying  the  consistency 
condition  at  a  finite  difference  point  using  a  backward  difference  formula  for  the  first 
derivative.  The  yield  function  differential  is  given  as: 


At  this  point  the 


'~'n+\ 


do„^,=0  (2.1-8) 

iteration  counter  k  is  introduced.  Values  at  A:  =  1  are  given  as: 

(2.1-9) 


(2.1-10) 


(2.1-11) 
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+  Yn^l  ^ 


trial 


Ak) 


M. 

da  ) 


\(k) 


n+1 


(2.1-12) 


For  elastic  steps,  the  yield  function  (Equation  2.1-10)  would  be  less  than  or  equal 
to  zero  and  the  trial  stress  would  be  substituted  as  the  new  stress.  If  the  yield  function 
were  greater  than  zero  the  Newton-Raphson  iteration  would  be  required  to  drive  the  yield 
function  (Equation  2.1-10)  and  residuals  (Equation  2.1-12)  to  zero.  It  is  interesting  to  note 
that  for  the  first  iteration  the  local  residuals  are  all  zero  except  for  the  yield  condition  (f) 
which  is  not  satisfied. 

A  Taylor  series  expansion  is  made  on  the  residuals  and  yield  condition  (Equations 
2. 1-12  and  2.1-10)  at  time  and  truncated  at  the  linear  term.  This  step  amounts  to  a 
derivation  of  the  Newton-Raphson  nonlinear  solution  algorithm.  Both  approximations  are 
set  to  zero: 


R. 


w 

n+1 


+  <‘J=0 


(2.1-13) 


To  obtain  accurate  results,  it  is  assumed  that  the  step  sizes  are  sufficiently  small  that 
the  differentials  can  be  approximated  with  finite  increments  (i.e.,  ==  and 

4fn+i  ==  This  is  required  to  assure  accuracy  of  the  finite  difference  approximations. 

In  addition,  care  must  be  taken  in  assuring  that  the  global  step  sizes  (and  the  step  sizes  used 
in  the  local  Newton-Raphson  iteration)  are  small  enough  so  that  there  is  convergence  at  the 
local  level.  This  is  often  accomplished  by  including  substepping  at  the  local  level. 

Equations  2.1-7  and  2.1-8  are  substituted  into  Equations  2.1-13  and  2.1-14, 
respectively,  and  Equation  2.1-13  is  multiplied  by  C ''  (assuming  it  exists).  The  resulting 
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equations  are  given  in  incremental  form  as: 


^n+I  ^  — «  +  l  ^^n+l  ^  ^/n  +  l  -  ^ 


(2.1-15) 


where  +ri^,\ 


f  V  ^ 

+  f-  A<,  =  0 

w<7y„+i 


(2.1-16) 


Note  that  the  H  matrix  is  shown  as  an  inverse  matrix  for  clarity  purposes  in  this  section. 

However,  in  the  coding  of  the  algorithm,  standard  matrix  factorization  and  substitution  is 
usually  used. 

Equation  2.1-15  can  be  rewritten  to  solve  for  the  unknown  increment  in  stresses 
given  an  increment  in  the  plasticity  parameter.  The  stress  increments  are  given  as: 


AC,  = 


“n+l  ^"■^‘^^'”■^'1^(7 


(2.1-17) 


The  increment  in  stresses  of  Equation  2. 1-17  can  be  substituted  into  Equation 


2.1-16  yielding: 


f(k)(  ^  Y  ^  =■(*)  +  Ay^*^  {  ^  1  -  0 

/n+l  ^«+l  -'Si+l^^/n+1 

\o(yJn^\ 


(2.1-18) 
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Finally,  Equation  2.1-18  can  be  rewritten  to  solve  for  the  increment  in  the  plasticity 
parameter  as: 


fW  _ 

Jn  +  \ 


Ar: 


(k)  _ 

+1 


da 


^n  +  \  ^ 


Vl  +  1 


Jn  +  \ 


yda 


T  (1. 

—n+l 

A+1  V 


do 


\(i) 


Jn+l 


(2.1-19) 


The  Closest  Point  method  begins  by  calculating  the  “trial”  elastic  stress  (Equations 
2.1-5  and  2.1-9).  This  is  tantamount  to  the  first  Newton-Raphson  iteration  starting  with  an 

initial  guess  of  an  elastic  step  (i.e.,  Acr^”^,  =  If  the  yield  function  is  greater  than 

zero,  this  stress  point  is  taken  as  the  first  iterate.  The  yield  function,  derivatives  and 
residuals  are  then  calculated  at  this  stress  point.  These  values  are  then  inserted  into 
Equation  (2. 1-19)  and  the  increment  in  the  plasticity  parameter  is  determined.  The 
plasticity  parameter  increment  is  inserted  into  Equation  2.1-17  to  solve  for  the  stress 
increments.  The  unknown  plasticity  parameter  and  stresses  are  then  updated  via: 


- 

/  n  +  l 


Yn^l+^Yn 


(k) 

+1 


(2.1-20) 


r(*+l)  _  ^(-k) 


r(k) 


'n  +  l 


=  CT;V,+Acr;,, 


(2.1-21) 


The  new  stress  point  and  plasticity  parameter  are  used  for  the  next  iteration  with  the 
process  continuing  until  the  yield  function  and  stress  residuals  are  within  some  tolerance  of 

zero.  Recall  (from  Equation  2.1-4)  that  the  plasticity  parameter  shown  here  is,  in 


fact.  At  y\„+i-  The  Closest  Point  stress  point  algorithm  is  described  in  Box  2.1-1  and  the 
behavior  during  iteration  is  shown  Figure  2.1-2. 


Box  2.1-1.  Fully  Implicit  Closest  Point  Stress  Point  Algorithm  (simple  model) 


1.  Given: 


£  £’’  C 

^n+1’ 


2.  Initialize: 


3.  Calculate  elastic  prediction:  cr'™' =  C  (e„+i  -  e;f ) 

4.  Calculate  yield  function  and  local  residuals: 


rf(k)  _  {k)  trial  ,  (k)  p 

^«+i  -  +7«+i 


5.  Test  for  convergence:  IF  (Xz+i  <  TOLERf  AND  <  TOLER^j  EXIT 

6.  Compute  the  derivatives  and  assemble  matrices: 


7.  Obtain  increment  in  plasticity  parameter: 


f  ^ 

A^(k)  ^  _ 

/i+l  /  ^  n(/:)  X  \(/:) 

I  ^  I  7:(^)  I  ^  1 


8.  Obtain  increment  in  stresses: 


= 


9.  Update:  riT  =  +  Af„%  =  <.  +  A(T<« 


10.  Set: 


k  ^k+  l,Goto  step  4 
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2.2  Derivation  of  Consistent  Tangent  Moduli  (simple  model) 


When  the  Newton-Raphson  algorithm  is  used  to  solve  the  global  nonlinear 
problem,  a  global  Jacobian  matrix  is  required.  This  global  Jacobian  matrix  is  an  assembly 
of  element  Jacobian  matrices  which,  in  turn,  are  composed  of  Jacobian  matrices 
constructed  at  the  quadrature  points.  The  accurate  construction  of  the  Jacobian  matrix  at  the 
quadrature  point  requires  a  matrix  of  derivatives  of  the  stress  component  relative  to  the 
assumed  strain  increments.  This  matrix  is  given  the  name  “consistent  tangent  stiffness 
matrix”  [Simo  and  Taylor,  1985]  and  must  be  supplied  by  the  material  routine.  Recall, 
however,  that  an  “exact”  Jacobian  is  not  necessarily  required  for  convergence,  although  it 
generally  improves  the  rate  of  convergence. 

Derivation  of  the  consistent  tangent  moduli  will  follow  a  framework  similar  to  the 
Closest  Point  stress  point  algorithm.  For  eventual  comparison  to  the  classical  continuum 
developed  moduli,  indicial  notation  is  introduced  in  this  section.  The  stress  equations 
(Equation  2.1-5)  are  given  in  indicial  notation  as: 


^  ^trial  .. 


ijkl 


A. 

J„+x 


(2.2-1) 


where 

As  with  the  stress  point  algorithm,  the  stress  equations  are  differentiated  at  time 


t  However,  since  the  tangent  moduli 


K  J„+, 


are  the  variation  in  stresses  given  a 


variation  in  strains,  the  strain  at  time  is  not  fixed.  The  differentiated  stress  equations 
can  now  be  given  as: 


7n+\ 


ijkl 


(2.2-2) 


^ijki 
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Both  sides  are  multiplied  by  C and  rearranged  to  solve  for  the  increment  in 
stresses.  The  stress  increments  are  given  in  terms  of  the  increment  in  strains  and  plasticity 
parameter  and  are: 


n+l 


^  1 

^+1  _ 

(2.2-3) 


where  = 


Qw  +  7, 


n+\ 


d^f 


dG^dG„ 


Jn+l. 


The  increment  in  stresses  is  inserted  into  the  discrete  consistency  condition 
(Equation  2.1-8).  The  resulting  expression  for  the  increment  in  the  plasticity  parameter 

(^7n+j)  given  in  terms  of  the  increment  in  strains.  This  is  given  as: 


= 


^  jf ' 

Jn+l 


H  dE 

mnop  (>p„^^ 


(  ^  ] 

(  ^  1 

^abcd 

n+l 

y^^cd  J 

n+l 

(2.2-4) 


Equation  2.2-4  is  substituted  into  Equation  2.2-3  which  now  puts  the  stress 
increments  in  terms  of  the  strain  increments: 
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""tiki 


dSui 


d£„„ 

tftnojy  + 1  f 


n+l 
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^ahcd 


A. 


K^^klJn^X 


(2.2-5) 


The  consistent  tangent  moduli  are  determined  by  differentiating  the  stress 


increments  by  the  strain  increments  at  time  The  moduli  are  given  as: 


^  do,  ^ 


ijmn 
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f  ^  1 

n+1 

'‘opkl 


^ijkl 


[  ^  1 

f  ^  ^ 

\do^^  j 

^ahcd 
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(2.2-6) 


where  = 


<^yr/+7n+l 


ydo^da„„y^^ 


-1 


The  tangent  moduli  for  the  associative  flow  rule  and  elastic-perfectly-plastic  material 
are  derived  using  the  classical  continuum  approach  [Dafalias,  1990]  independent  of  any 
particular  stress  point  algorithm.  The  resulting  moduli  are  given  as: 
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D  =C 

^ijkl  ^ijkl 
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df  1 

\.do„„  J 

J 

’'opkl 


r  ^ ' 

{  ^  ] 

^  abed 

(2.2-7) 


where  D-^j^  -  tangent  moduli  derived  from  continuum  plasticity. 

The  consistent  tangent  moduli  developed  by  the  Closest  Point  algorithm  (Equation 
2.2-6)  results  in  terms  that  are  consistent  with  the  method  used  to  integrate  the  stress  point 
algorithm.  Upon  comparing  the  moduli  it  is  seen  that  the  difference  between  the  two 

approaches  for  this  simple  model  lies  in  the  use  of  the  E  matrix  instead  of  the  elastic 

constants.  The  H  matrix  depends  not  only  upon  the  elastic  moduli  but  also  the  plasticity 

parameter  and  the  second  derivative  terms  of  the  yield  function.  Earlier  literature  has 
shown  that  consistent  tangent  moduli  usually  result  in  faster  global  convergence  [Simo  and 
Taylor,  1985]  because  they  create  an  exact  global  Jacobian  (unless  other  approximations 
are  made).  Reduction  in  global  calculations  can  result  in  significant  cost  savings  if  the 
increased  local  calculations  are  not  overwhelming. 
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2.3  Derivation  of  Closest  Point  Algorithm  (general  model) 

Classical  plasticity  includes  internal  variables  that  describe  the  changes  in  shape, 
size  and  location  of  the  yield  surface  as  plastic  straining  occurs.  The  general  formulation  of 
the  Closest  Point  algorithm  was  developed  with  the  internal  variables  expressed  in  rate 
form  [Simo  and  Hughes,  1990].  The  yield  function  for  the  general  formulation  is  now 
described  in  terms  of  the  stresses  and  internal  variables  and  is  given  as: 

/  =  /{ff,?)  =  0  (2.3-1) 

where  q  =  internal  variables. 

The  rate  equations  for  the  plastic  strains  with  non-associative  flow  and  internal 
variables  are  given  in  a  genercd  form  as: 

£'’=YN  (2.3-2) 

q  =  Y  h  (2.3-3) 

where  N  =  N(o,  q)  =  direction  of  plastic  strains 

h  =  h(G,  q)  =  direction  of  internal  variables. 

The  plastic  strain  and  internal  variable  directions  are  shown  as  functions  of  both  stresses 
and  internal  variables,  although  they  can  be  a  function  of  just  one  of  these  characteristics. 


The  rate  equations  for  both  the  strains  and  internal  variables  are  approximated 


between  times  and  using  a  backwards  Euler  scheme  and  are  given  as; 

(2.3-4) 

‘S'n+I  ~  “*“Tn+l^n+l  (2.3-5) 

where  =  N( q^^,) 

^n  +  l  ~  M^n  +  1<  ^n+l)- 


As  in  the  stress  point  algorithm  (Section  2.1),  the  plastic  strains  at  time  )  are 

added  and  subtracted  in  the  stress  equation  (Equation  2.1-3)  and  use  is  made  of  Equation 
2.3-4.  The  resulting  equation  for  stress  is  in  terms  of  the  “trial”  stress  (elastic  prediction) 
and  is  given  as: 

(2.3-6) 


where  a'™' =  C(e„„ -e^). 


The  set  of  nonlinear  equations  to  be  solved  consists  of  the  stress  equations 
(Equation  2.3-6)  and  the  internal  variable  equations  (Equation  2.3-5).  The  local  residual 
vector  for  these  equations  is  defined  as: 


-Qn  -Yn^A^^ 


(2.3-7) 


The  local  residuals  are  differentiated  at  time  to  give  the  local  Jacobian  matrix. 
Recall  that  the  strain  rates  and  internal  variable  rates  can  be  functions  of  both  stress  and 
internal  variables.  Again,  since  the  strains  are  given  for  the  stress  point  algorithm  they  are 


fixed  at  time  t^^j.  The  differentiated  residuals  are  given  as: 


d(yn.i+dr„^,CN„,,+y„,,C 


V<?CTj„^, 


^^n+l  n+\^n+l  ^*^n+l  Y  n+\ 


\dqj 


\dqj 


dq 


w+1 


«+l 


n+] 


=  0  (2.3-8) 


The  yield  function  is  differentiated  at  time  and  now  includes  the  contributions  of 
both  stress  and  internal  variables.  The  discrete  consistency  condition  is  now  given  as: 
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dC^n^l  + 


[dqj 


dqn^x  =0 


n+l 


(2.3-9) 


To  facilitate  the  general  algorithm  development,  a  number  of  matrices  and  vectors 
are  defined: 


G  = 
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0 
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-I 
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’dN 
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'1  0 
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da 

dh 

dq 

dh 
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0  1  ••• 

.dq. 

da 

dq. 

- 

(2.3-10) 


Equations  2.3-8  and  2.3-9  can  now  be  written  in  matrix  form  as: 

a  Z«>  +y»>  G  <«,  =  0  (2.3-1 1) 

V/“;  =  0  (2.3-12) 

Substituting  Equations  2.3-1 1  and  2.3-12  into  the  linear  portion  of  the  Taylor  series 


expansion  and  assuming  a  finite  increment  for  the  differentials  (Equations  2.1-15  and 
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2.1-16)  results  in: 

+  AZ™,  +  Ay»  G  Z®  +r"’,  G  Z®  AZ®,  =0  (2.3-13) 

/®-hV/®AZ®=0  (2.3-14) 

Equation  2.3-13  can  be  solved  for  the  increment  in  stresses  and  internal  variables 
yielding: 

[G-C,  +  ]  (2.3-15) 

where  E  =  [G"‘+rf+i  y„+|]  . 

As  noted  in  Section  2.1,  the  inverse  of  the  E  matrix  is  shown  for  clarity,  however, 

standard  matrix  factorization  and  substitution  would  generally  be  used  in  the  algorithm. 
Finally,  the  increment  in  the  plasticity  parameter  is  determined  by  substituting  Equation 
2.3-15  into  Equation  2.3-14  and  is  given  as: 


= 


f(k)  _  ^Ak)  2  ^-1  J^{k) 

Jn+l  ^Jn  +  l  ^  +  \ 

Vn-H  “  ^n+\ 


(2.3-16) 


As  before,  the  Closest  Point  algorithm  begins  with  the  calculation  of  the  “trial” 
stresses  (Equation  2.3-6).  If  the  yield  condition  is  positive,  the  derivatives  and  residuals 
are  calculated  at  this  stress  point.  The  increment  in  the  plasticity  parameter  (Equation  2.3- 
16)  is  evaluated,  then  the  increment  in  stresses  and  internal  variables  are  calculated  from 
Equation  2.3-15.  The  plasticity  parameter,  stresses  and  internal  variables  are  then  updated 
via: 


/ n-f-1  /  n+\  *  n+1 


y(*+l)  _  y(k)  ,  Ay(k) 
^n+l  ~  ^n+1  ^ 


(2.3-17) 


(2.3-18) 
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The  new  stress  point,  internal  variables  and  plasticity  parameter  are  used  for  the 
next  iteration  with  the  process  continuing  until  the  yield  function  and  the  local  residuals  are 
within  some  tolerance  of  zero.  The  general  stress  point  algorithm  is  described  in  Box  2.3-1 
and  shown  graphically  in  Figure  2.3-1. 

Box  2.3-1.  Fully  Implicit  Closest  Point  Algorithm  with  Internal  Variables 


k  =  l,  rl'i,=o,  qZ=Qn 


1.  Given:  e„+i,  e^,  C 

2.  Initialize: 

3.  Calculate  elastic  prediction:  cr'""' =  C  (e„+i  -  ) 

4.  Calculate  yield  function  and  local  residuals: 


^n+l  “ 


.^n+l  ~  qn~  7n+l^n+l 


5.  Test  for  convergence:  IF  <  TOLERf  AND  <  TOLER^j 

6.  Compute  derivatives  and  assemble  matrices: 

2l*.’,=[G-'+y“  !'®r 

7.  Obtain  increment  in  plasticity  parameter: 


EXIT 


A7r.i 


f  _  Vf 

(k)  _  Jn+\  ^Jn+\  ^n+\  ^  +  I 


-::U)  7(^) 

Vn+1  ^n+\  ^n+\ 


8.  Obtain  increment  in  stresses  and  internal  variables: 

A£l‘.>.=-E™  [G-'C  +  Ar“>,Z;‘.’,] 

9.  Update;  r“;"  = /“  +  AT®,.  E" V  =  S®,  +  AX“> 

10.  Set:  A:  ^  +  7,  Go  to  step  4 
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2.4  Derivation  of  Consistent  Tangent  Moduli  (general  model) 


Derivation  of  the  consistent  tangent  moduli  is  somewhat  more  complicated  because 
the  internal  variables  are  now  included.  The  stress  and  internal  variable  approximations 
(Equations  2.3-6  and  2.3-5)  are  given  as: 


(2.4-1) 


«..i  (2.4-2) 

where  <7"=C(£„,-£r). 

Equations  2.4-1  and  2.4-2  are  differentiated  at  time  Recall  that  the  tangent 
moduli  are  variations  in  stress  with  variations  in  strain  and  therefore  the  strains  at  time  ^n+I 
are  not  fixed.  The  stress  and  internal  variable  differentials  are  now  written  as: 


n+\  n+\  ^  ^n+\  Y n+t  ^ 


JdN]  , 

■Yn^iC  —  dq„,^ 

V  ^  )n+\ 


(2.4-3) 


,  ,  (dh\  \dh\  n  A  A\ 

Solving  Equations  2.4-3  and  2.4-4  for  the  increments  in  stresses  and  internal 
variables  and  writing  in  matrix  form  results  in: 

■E..,  [dJC  +  drZL,  (2.4-5) 


where  X  = 


Equation  2.4-5  can  be  substituted  into  the  discrete  consistency  condition  (Equation 


2.3-12)  and  solved  for  the  increment  in  the  plasticity  parameter  as: 


v/„.. 


^n+X  -^n+l 

E  7 

—n+X  ^n+X 


(2.4-6) 


Substituting  the  plasticity  parameter  (Equation  2.4-6)  into  Equation  2.4-5  results  in 
an  expression  for  the  increments  in  stress  and  internal  variables  in  terms  of  strains. 
Differentiating  this  expression  with  respect  to  strains  at  time  ^n+]  provides  the 
algorithmically  consistent  tangent  moduli  for  the  general  model: 
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where  = 


C-^+Yn.X 

Yn.X 


Yn.X 

ydq) 

(  dh\ 
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-i+r, 
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Jn+\ 


-1 


Since  the  tangent  moduli  describe  how  the  stress  increment  changes  in  relation  to 
the  strain  increment,  the  5  matrix  is  partitioned  (because  it  contains  terms  for  both  the 
stresses  and  internal  variables).  It  is  important  to  note  that  the  tangent  moduli  are  generally 
not  symmetric  since  5,2  and  Ej,  are  generally  not  equal.  Because  of  the  general  lack  of 

symmetry,  the  entire  matrix  needs  to  be  available  and  nonsymmetric  solution  methods  are 
required  at  the  stress  point  algorithm  level.  It  is  worth  noting  that  most  “real  world” 
problems  (e.g.,  large-deformations,  slide-surfaces,  etc.)  are  nonsymmetric  even  if  the 
constitutive  relations  are  symmetric. 

The  unique  case  for  symmetry  comes  when  considering  associative  flow  (for 
strains)  and  what  is  coined  as  “associative  hardening”  (for  internal  variables).  These  are 


given  respectively  as: 


^  =  1^  h  =  ^  (2.4-8) 

3a  dq 

It  must  be  noted  that  this  symmetry  is  only  possible  when  the  elastic  behavior  is  derivable 
from  an  elastic  potential  so  that  the  C  matrix  is  symmetric.  This  condition,  in  conjunction 

with  Equation  2.4-8,  results  in  a  symmetric  E  matrix,  but  introduces  second  derivatives  of 
both  the  yield  function  and  the  internal  variable  rate  equations.  The  symmetric  E  matrix  is 


given  as: 


(2.4-9) 


Substituting  this  matrix  into  the  tangent  moduli  (Equation  2.4-7)  results  in 
symmetric  tangent  moduli.  Most  material  models,  however,  rarely  include  this  “associative 
hardening”  and  often,  such  as  the  case  of  the  soil  model  considered  here,  do  not  have  a 
symmetric  elastic  matrix  (C). 

The  tangent  moduli  for  the  associative  flow  rule  and  elastic -perfectly  -plastic 
material  using  the  classical  continuum  approach  [Dafalias,  1990]  is  given  as: 
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Kp  + 


\daj 


da  ) 


(2.4-10) 


where 


ydq) 


h. 


The  differences  between  Equations  2.4-7  and  2.4-10  for  the  more  general  model  is 
not  quite  as  simple  as  was  the  case  in  Section  2.2.  The  algorithmically  consistent  moduli 
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include  additional  derivatives  from  both  the  strain  and  internal  variable  rate  equations.  In 
the  continuum  approach,  however,  the  internal  variables  enter  into  the  moduli  only  through 
the  hardening  modulus  {K^  in  the  denominator;  whereas  in  the  consistent  moduli,  they 

appear  in  both  the  numerator  and  the  denominator  through  the  S  matrix  and  the  Z  vector. 

An  approach  for  creating  symmetric  tangent  moduli  is  given  in  the  following 

section. 

2.5  Derivation  of  Closest  Point  Algorithm  (implicit-explicit  model) 

An  alternate  formulation  of  the  Closest  Point  algorithm  is  developed  similar  to  the 
general  model,  but  involves  treating  the  internal  variable  rates  in  an  explicit  fashion.  The 
rate  equations  for  the  plastic  strains  are  approximated  between  times  to  K+i  as  before 
using  a  backward  Euler  form: 

(2.5-1) 

where 

The  internal  variable  rates,  however,  are  approximated  with  the  plasticity  parameter 
evaluated  at  time  but  the  direction  evaluated  at  time  4.  This  is  given  as: 

(2.5-2) 

where  h„  =  h( cr„,  qj. 

The  local  residual  vector  is  defined  as: 

Jln+\  y  n+\^n 


(2.5-3) 
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The  local  residual  is  differentiated  at  time  Recall  that  the  internal  variable  rates 
are  evaluated  at  time  4  and,  therefore,  have  no  contribution.  The  differentiated  residuals 
are  given  as: 


dK.,  = 


dq„^i  -dr„^A 


\d<y)n^l 


^CS’n+I  +7n+l<^ 


V  ^‘1  Jn^\ 


dq 


n+l 


=  0  (2.5-4) 


Equation  2.5-4  is  solved  for  the  increment  in  stresses  and  internal  variables 
(Equation  2.3-15)  and  substituted  into  the  discrete  consistency  condition  (Equation  2.3-9). 
The  increment  in  stresses,  internal  variables  and  the  plasticity  parameter  are  given  as: 
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where  E  =  [G  ‘  -i-  y  t] 
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0 

0 


(2.5-6) 


The  Y  matrix  with  all  the  zeros  is  shown  for  discussion  purposes.  In  the  actual  algorithm, 
only  the  nonzero  submatix  would  be  stored  and  the  number  of  actual  calculations  would  be 
reduced  accordingly. 

The  consistent  tangent  moduli  are  obtained  as  in  Section  2.4.  Differentiating  this 


expression  with  respect  to  strains  provides  the  consistent  tangent  moduli: 
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where  E  = 


C*'  +r 


n+\ 


V<?CTj„+, 


n-i 


Symmetric  tangent  moduli  are  obtained  when  both  the  strain  rate  functions  and  the 


r 


elastic  moduli  are  symmetric 


i.e.,  associative  flow,  N  =  — 

da  j 


.  Because  the  internal 


variables  are  evaluated  at  time  “associative  hardening”  is  not  required  for  symmetry. 
Comparing  these  moduli  to  the  continuum  moduli  provides  an  interesting  insight.  The 
continuum  moduli  for  non-associative  flow  are  defined  as: 


CN 


D  =  C- 


\daj 


da 


C  N  +  K„ 


(2.5-8) 


where 


h. 


yoqj 


The  consistent  moduli  appear  similar  to  the  continuum  derived  moduli  (Equation 
2.5-8)  except  for  the  second  derivative  terms,  which  are  included  with  the  elastic  terms  of 
the  consistent  moduli.  The  contribution  of  the  internal  variables  now  appears  in  the 
denominator  for  both  methods.  The  difference  is  that  the  hardening  modulus  (K^)  is 
defined  at  time  for  the  continuum  moduli  and  at  time  for  the  implicit-explicit 
consistent  moduli. 

An  additional  modification  can  be  made  by  updating  the  directions  that  are  evaluated 


at  time  by  directions  calculated  using  the  previous  iterate  as: 

NL=N(o'::ur)  p-s-?) 

P'5-IO) 

2.6  Application  of  the  Closest  Point  Method  for  Bounding  Surface 
Plasticity 

The  Closest  Point  method  was  developed  in  the  previous  sections  for  typical 
plasticity  models  that  are  elastic  within  the  yield  surface.  The  Bounding  Surface  Plasticity 
model,  however,  uses  the  surface  not  as  a  yield  surface,  but  as  a  bound.  The  complication 
in  Bounding  Surface  Plasticity  is  that  yielding  can  occur  within  the  bound.  With  a  clever 
choice  of  internal  variables  the  model  can  be  set  up  to  resemble  a  classical  yield  model. 

This  section  will  describe  the  establishment  of  the  functions  and  their  derivatives  for  the 
Bounding  Surface  Clay  model. 

The  approach  used  for  this  development  is  to  separate  the  volumetric  and  deviatoric 
components  of  the  stress  and  strain  relationships  similar  to  the  Cam-Clay  model  [Meschke 
and  Simo,  1994].  This  allows  for  developing  the  relationship  between  the  volumetric 
stress  (7 )  and  the  Bounding  Surface  size  (/{,),  which  will  reduce  the  number  of  nonlinear 
equations.  The  nonlinear  elastic  volumetric  relationship  is  described  in  Section  1.3. 
Because  of  the  Macaulay  brackets  in  Equation  1.3-1,  there  are  two  elastic  regions 
delineated  by  the  transition  volumetric  stress  (/;).  In  order  to  remain  consistent  with 
Section  1.3,  the  equations  in  this  section  where  /  >  7,  are  again  denoted  as  the  (a)  equations, 
and  the  equations  where  7  <  7,  are  again  denoted  as  the  (b)  equations.  For  this  section,  it  is 
assumed  that  the  steps  that  cross  the  transition  volumetric  stress  are  sufficiently  small  that 
the  elastic  relationships  for  each  region  will  provide  an  adequate  approximation  for  the  step. 
Separate  equations  for  steps  crossing  the  transition  stress  are  developed  later  in  Section 
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3.2.  Rewriting  the  differential  relationships  of  the  volumetric  strain  and  stress  (Equation 
1.3-7)  in  terms  of  the  total  and  plastic  strains  yields; 

=  (2.6-ia) 

=  (2.6- lb) 

As  described  in  the  previous  sections,  to  formulate  the  initial  “trial”  stress  to  start 
the  Newton-Raphson  iteration,  an  elastic  prediction  is  made  assuming  no  plasticity.  During 
the  plastic  correction,  the  total  volumetric  strain  is  held  constant.  This  results  in  the  change 
in  volumetric  stress,  during  the  plastic  correction,  being  written  in  terms  of  only  the  plastic 
volumetric  strain.  This  is  given  as: 


^n+l 


where 


del, 

^n+\  P  ^n+l 


{l>h) 

(/</,). 


(2.6-2) 


Because  the  total  volumetric  strain  increment  is  fixed  and  applied  in  the  first  iteration,  the 
correction  which  involves  the  developing  plastic  volumetric  strain  is  negative.  It  must  be 
noted  that  during  the  plastic  iteration  some  plastic  volumetric  strain  will  be  created  and  thus 
change  the  elastic  prediction  in  the  next  iteration  (see  Box  2.6-1). 

The  relationship  for  the  internal  variable  {Q  of  the  Bounding  Surface  model  is 
developed  in  Section  1.5.  The  differential  relationship  is  given  by  Equation  1.5-7  and  can 
be  written  at  time  as; 

= «...  (Q 


where 
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Note  that  during  the  plastic  iteration  both  the  volumetric  stress  and  Bounding 
Surface  size  are  related  to  the  plastic  volumetric  strain  (Equations  2.6-2  and  2.6-3). 
Therefore,  the  differential  in  the  Bounding  Surface  size  can  be  written  in  terms  of  the 
differential  in  the  volumetric  stress  during  the  plastic  iteration  as: 

dl^  (2.6-4) 

The  plastic  strain  rates  for  the  Bounding  Surface  model  are  defined  in  terms  of  the 
plasticity  parameter  and  directions  given  by  the  “image”  stresses  (Equation  1.2-2). 
Rewriting  these  rates  in  terms  of  volumetric  and  deviatoric  strains  gives: 


(2.6-5) 

II 

(2.6-6) 

where 

plastic  deviatoric  strain  tensor 

s  =  bs  =  deviatoric  “image”  stress  tensor. 

The  Bounding  Surface  function  (F)  used  for  this  study  is  the  single  ellipse  model 
[Kaliakin,  1985;  Herrmann,  et.  al.,  1985;  Herrmann,  et.  al.,  1987]  as  described  in  Section 
1.4  (Equation  1.4-1). 

As  in  the  previous  sections,  the  increments  are  approximated  between  times  and 
with  the  backward  Euler  approximation,  which  results  in; 


-  QP 


dl 


J  n-^\ 


(2.6-7) 


=<+Yn.^ 


(2.6-8) 
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Recall  that  the  notation  y„^j  will  be  used  instead  of  the  more  appropriate  At  for  clarity. 

The  deviatoric  stress  rates  are  defined  in  terms  of  the  deviatoric  strain  rates  and  the 
shear  modulus  (G)  at  time  and  are  given  as: 

s  =-2G  e’’ ,  (2.6-9) 

The  rate  equation  for  the  loading  surface  (b)  was  developed  in  Section  1.7.  For  this 
formulation  of  the  Closest  Point  algorithm,  the  variable  b  will  be  treated  as  an  additional 
internal  variable.  From  Equation  1.7-1 1,  a  backward  Euler  approximation  yields: 

=K+rn^l8n^\  (2.6-10) 


where 


8n+l  ~ 


K,-bK^-3Cv„  {\-b) 


L  (dF 


\dl  j 


X-K 


^2  A 


,  .dF  dF 


'  n+\ 


The  local  residuals  for  the  nonlinear  equations  are  defined  using  the  approximations 
of  the  plastic  deviatoric  strains  (e),  the  plastic  volumetric  strain  (ff)  and  the  loading  surface 
(b).  Rewriting  Equations  2.6-7,  2.6-8  and  2.6-10  gives: 


^dF^ 


\dsj 


(2.6-1  la) 


n  +  l 


dF 

dJ 


n+l 


(2.6-1  lb) 


X+l  =^n+I-^n-rn+1  Sn+I 


(2.6-1 Ic) 


where  R  =  residual  equations  for  the  six  deviatoric  strain  expressions 
=  residual  equation  for  the  volumetric  strain  expression 
-  residual  equation  for  the  loading  surface  (b)  expression. 
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Since  the  internal  variable  for  the  Bounding  Surface  size  {1^)  is  expressed  as  a  function  of 
the  plastic  volumetric  strain  (Equation  2.6-3),  it  will  be  incorporated  later  in  this  section  and 


will  not  be  treated  as  a  separate  residual. 

As  in  the  general  Closest  Point  framework,  the  residuals  are  differentiated  at  time 
noting  that  E  is  a  function  of  the  stresses  /  and  s  and  the  internal  variables  4  and  b. 
The  differentials  of  the  residuals  are  given  as: 


dsds 


^^n+l  ^n+\  Tn+I 


ds 


n+l 


■Yn 


+1 


-r..i 


dsdb 


-  dYn^x 


dl. 


n+l 


■Yn.l 


^  d^F  ^ 


n+l 


dsdi 


dr 


«  Jn+\ 


J  n+l 


(2.6- 12a) 


= 


del,  -Yn.X 

^  d"-F^ 


Y 


n+] 


didb 


^  d^F^ 

db,,,  -  dr,,, 


J  n+l 


d^n+l  Y  n+l 
V  dl 


d^F 

dldl 


dln+l  -Yn+l 


^  d^F  ^ 


J  n+l 


dldl 


dl. 


o  Jn  +  \ 


(2.6-12b) 


'dR,„  =  db„„  -  r„„ 

'"dg^ 


dS„„  -  Yn+l 


-Yn+l 


\dbj,„ 


\d^)n+l 

db,„  -  dr  rt+1  Sn+] 


I'  ^ 
\dl  j,„ 


dr+\  Y  n+l 


'  dg^ 


dr 


(2.6-12C) 


The  differentials  of  the  plastic  volumetric  and  deviatoric  strains  (i.e.,  del,,ddl,) 
can  be  expressed  in  terms  of  stresses  by  substituting  Equation  2.6-2  and  the  differential 
form  of  Equation  2.6-9  into  Equation  2.6-12.  The  internal  variable  4  is  incorporated  by 
substituting  Equation  2.6-4  into  Equation  2.6-12.  The  substitutions  result  in  equations  that 
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are  expressed  entirely  in  terms  of  stresses  and  internal  variables  and  are  given  as: 


dR„^x  = 


2  r?  ^ 


2G 


+r 


d^F 

^ds 


r  a2 


ds 


n+l  T«+l 


Jn+\ 


d'^F  H  d'^F 


^dl  K  ^dl, 


dl. 


n+l 


"  /n  +  l 


^  1)2  17  \ 


■r„+i 


d^F 

dsdb 


db. 


/l  +  l 


dy 


rt+1 


J  n  +  \ 


/n+I 


(2.6- 13a) 


'dK,,  =  -  r„+i 


-r 


n+l 


dids 
d^F 


^n+l  -  T7  +  r 


V^^^'^'^yn+l 

/  n2  rr  'N 


// 

ydidl  K  dl dig  j 


dl. 


n+I 


n+l 


dldb 


dbg,,  -dy. 


y  n+l 


n+I 


dl  y„+i 


(2.6- 13b) 


dRg„  =-r„+i 
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\dSjg„ 


d^n+\  y  n+] 


ydl  Kdlgj 


dl. 


n+l 


n+l 


+ 


dg^ 


db. 


n+I 


dy„„  gn+l 


/n+l 


(2.6- 13c) 


The  Bounding  Surface  function  (F)  is  differentiated  noting  that  it  is  a  function  of  s, 
1,  /q  and  b.  This  differentiation  results  in  the  discrete  consistency  condition  which  is  given 


as: 


dF  = 


ds  Jn^l 


ds„„  + 


V  dl  y„„ 


-I 


dl  + 


V^^^n^n+l 


Kdb 


/n+I 


A+1  =  0 


(2.6-14) 


Recalling  the  relationship  between  7  and  Ig  (Equation  2.6-4).  Equation  2.6-14  can 
be  rewritten  as: 


dF  = 


V  ds  Jg„ 


ds_^,  + 


^dF  H  dF^ 


dl  K  dl. 


o  Jn+\ 


(dF'^ 
dln^\  +  3- 


dbg,,  =  0 


(2.6-15) 
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The  residual  differentials  and  discrete  consistency  condition  can  be  written  in  matrix 
form  as: 


^^n+\  ~  ^n+\  ^^n+\  ^7 n+l 


n+1 


=  0 


where  =  [s  7  Bf^_ 

'dF  dF 


n  +  1 


^n+I  ■“ 


^  dl  ^ 
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= 

''^n  +  l 


dF 


dF  H  dF 


dF 
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ds 

~KdlJ 

db 

J/i+i 

(d^F^ 

f  d^F 

H  d^F  ^ 

J 

7 
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K^dlJ 
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^  d^F^ 
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r  32 
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d^F  H  d^F 
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dIdJ  Kdidl, 
dg  H  dg^ 
dl  Kdl„j 


oj 


(2.6-16) 

(2.6-17) 


f  32  17  A 


-7 

-7 

1-7 


d^F 
Sdb 
^  d^F^ 


ydidb  j 
\dh) 


In+I 


The  Jacobian  matrix  (A)  can  be  simplified  by  noting  that  the  second  derivatives 
involving  orthogonal  directions  (i.e.,  7  and  s  or  7^  and  s)  are  zero.  The  matrix  is  reduced 
to: 
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^ds 
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2G 
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«+i  (2.6-18) 
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The  second  derivatives  of  the  Bounding  Surface  function  (F)  are  derived  in 
Appendix  A.  Because  of  the  complex  form  of  the  b  rate  equation  (Equation  2.6-10),  the 
derivatives  of  g  are  approximated  with  finite  difference  formulas. 

As  in  the  previous  sections,  both  the  residuals  and  consistency  condition  are 
linearized  (i.e.,  Taylor  series  expansion  to  the  linear  term)  and  set  to  zero.  The  assumption 
is  made  that  step  sizes  are  sufficiently  small  so  that  the  differentials  can  be  approximated 
with  finite  increments.  The  resulting  equations,  including  an  iteration  counter  (k),  are 
given  as: 


+ {a-']Z  z“,  =  0 


(2.6-19) 


C  +  VF<»a4"=0  (2.6-20) 

Equation  2.6-19  is  rearranged  to  solve  for  the  unknown  increments  in  s,  I  and  b, 
and  is  given  as: 

AX*.’,  =  -A™,  [«+  Ar  Z]“  (2.6-21) 

Equation  2.6-21  is  substituted  into  Equation  2.6-20  and  solved  for  the  increment  in 
the  plasticity  parameter  (Af),  resulting  in: 


A  y(^)  n+\ 

^in+\  X 


(2.6-22) 


Once  the  increment  in  the  plasticity  parameter  is  obtained,  the  increments  in  the 
stresses  (s  and  7)  and  internal  variables  (b)  are  found  with  Equation  2.6-21.  The  internal 
variable  (Ig)  representing  the  Bounding  Surface  size  is  indirectly  included  via  the  plastic 
volumetric  strain  (Equation  2.6-3).  The  increments  in  plastic  volumetric  and  deviatoric 
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strains  are  calculated  using  Equations  2.6-2  and  2.6-9  and  are  given  as: 


1 

< 

1 _ 

{^) 

1 

1  ’ 

W 

As 

_Ar_ 

n  +  1 

_2G 

~K 

n+\ 

AI 

r 


J;i+1 


(2.6-23) 


where  1  =  [l  1  1  1  1  l]  • 

The  plastic  strains  and  internal  variable  {b)  are  then  updated  during  the  iteration  by: 

(2.6-24a) 


r"*'' 

pP  =  4-  Ap^ 


^n  +  \ 


=  e::;+AeC 


(2.6-24b) 


Jk) 


^n+1  ^n  +  \  ^  ^n+! 


(2.6-24C) 


The  elastic  “trial”  stress  is  calculated  with  the  plastic  strains  initially  set  to  the  values 
at  the  end  of  the  last  global  step.  The  values  of  I  and  I„  are  calculated  with  Equations  1.3-8 
and  1.5-8,  respectively.  The  deviatoric  stress  (s)  is  calculated  at  time  by: 


s 


n+] 


(2.6-25) 


The  Bounding  Surface  function  (F)  is  evaluated  with  the  “trial”  stresses  and  the  internal 
variable  (b)  initialized  with  the  value  at  time  If  F  is  less  than  or  equal  to  zero,  the  step  is 
elastic  and  the  new  stresses  are  sent  back  to  the  main  routine.  It  is  important  to  note  that  the 
internal  variable  (b)  is  a  measure  from  the  current  stress  point  to  the  Bounding  Surface  and 
has  to  be  recalculated  elastically  (not  evolved),  even  for  an  elastic  step.  If  F  is  greater  than 
zero,  the  trial  stress  state  is  used  as  the  initial  guess  for  the  plastic  iteration.  A  case  that  is 
unique  to  the  Bounding  Surface  model  occurs  when  the  elastic  “trial”  step  is  outside  of  the 
elastic  bound  but  does  not  go  outside  of  the  Bounding  Surface.  Treating  the  stress  measure 
(b)  as  an  internal  variable  allows  for  proper  evaluation  of  the  consistency  condition.  Using 

the  “trial”  stresses  and  b  evaluated  at  time  results  in  the  image  stress  state  )  being 
outside  of  the  bound  (F  >  0),  thus  indicating  plasticity.  This  is  shown  graphically  in 


Figure  2.6-1.  The  plastic  iteration  begins  evolving  the  plastic  strains  and  the  internal 
variable  (now  using  Equation  2.6-10)  and  continues  until  the  value  F  is  less  than  some 
tolerance  associated  with  the  Bounding  Surface  function  and  the  norm  of  the  residuals  is 
less  than  some  tolerance  associated  with  the  residuals.  At  convergence  the  image  stress 
state  coincides  with  the  bounding  surface.  The  plastic  algorithm  is  shown  in  Box  2.6-1. 
The  algorithmically  consistent  tangent  moduli  are  the  same  as  described  in  Section  2.4 
(Equation  2.4-7). 
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Box  2.6-1 .  Fully  Implicit  Closest  Point  Algorithm  for  Bounding  Surface  Plasticity 


1.  Given: 

2.  Initialize:  k  =  0,  ri!’,  =  0,  C  =  €!  =  < 

3.  Evaluate  /,  I„,  s  and  associated  terms: 


=  3^  4. 


*^n+l  ” 


2G..,  ( 


^n+l  ^/i+l  , 


4.  Evaluate  the  functional  and  residuals: 


Tik)  j(k)  i(k)\  ^  QP  ^QP  7^^) 


K\-K 


5.  Check  convergence:  If  (  tolerance  p  and  (  tolerance,^  Then:  EXIT 

6.  Compute  derivatives  and  assemble  matrices  (note,  derivatives  of  g  are 
estimated  via  finite  difference): 


£L 

2G 


(  B^F^ 

0 

'd^F  H 

tC  b]  BI  ^^  ^ 

-r  — 

[b1^^ 

\-y\-L\ 

labj 

dl  K  dl„) 

7.  Solve  for  increment  in  stresses  and  internal  variables, 

Z7(^)  _  Y7T?(^)  A(k)  Tf(k) 

\^Ak)  _  ^n+1  Az+l^*Vhl  A  r  I?  a-  An/ 


\7¥7(^)  A(k)  rw{k) 
^■^n+1  ^n  +  1  ^n+1 


=-/!»[/!  + Ay  Z]; 


8.  Solve  for  plastic  strains: 


1  J_ 
2  G  K 


9.  Update:  eC"  =  e„C  +  A^„C ;  0„C'’  =  =  *‘1!  +  A*<t> ;  r't;"  =  rf.’,  +  Af::. 

10.  Set:  k  <—  k+1  and  go  to  step  3. 


Bounding  Surface 
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3.  Reduced  Newton  Method 

The  Reduced  Newton  method  was  developed  by  Leonard  Herrmann  for  integrating 
the  two-dimensional  Cam-Clay  model  [Roscoe  and  Borland,  1968].  The  method  maps  the 
stress  rate  equations  and  the  consistency  condition  into  a  single  nonlinear  equation  and 
integrates  them  with  backwards  Euler  integration  using  Newton-Raphson  iteration.  The 
method  incorporates  uniform  substepping  and  has  been  implemented  into  a  steady-state 
finite  element  program  [Herrmann,  1997]. 

For  simplicity,  the  application  of  the  method  will  be  restricted  to  the  case  where 
there  is  no  dependency  on  the  “lode”  angle  (it  is  assumed  that  the  results  of  compression 
and  extension  triaxial  tests  are  approximately  the  same).  This  special  case  will  be  referred 
to  as  a  two-invariant  model. 

Applying  the  Reduced  Newton  method  to  the  two-invariant  Bounding  Surface 
Plasticity  model  for  clays  requires  a  second  equation  that  determines  the  level  of  plasticity 
occurring  within  the  Bounding  Surface  (Equation  1.2-4).  This  section  describes  the 
method  applied  to  the  two-invariant  version  of  the  clay  model. 

Section  3. 1  simplifies  the  Bounding  Surface  model  for  clays  to  a  two-invariant 
form.  Section  3.2  describes  elastic  stress  paths  that  begin  inside  the  elastic  nucleus  and 
intersect  with  it.  The  elastic  algorithmic  consistent  tangent  moduli  are  developed  in  Section 
3.3.  Section  3.4  explains  the  “Reduced  Newton”  stress  point  algorithm  for  plastic 
calculations.  The  corresponding  plastic  algorithmic  consistent  tangent  moduh  are 
developed  in  Section  3.5.  Finally,  Section  3.6  discusses  the  handling  of  the  unloading 
stress  path. 


3.1  Two  Invariant-Bounding  Surface  Plasticity 


Section  1  describes  the  Bounding  Surface  Plasticity  model  for  clays  when  all  three 
stress  invariants  are  involved.  The  model  can  be  simplified  when  it  is  only  a  function  of 
two  invariants.  The  Bounding  Surface  becomes  in  this  case: 

F  =  ^r+il-pl\"--ll=Q  (3.1-1) 

a, 

where  F  -  Bounding  Surface  function 
I=b(I-IJ+l, 

I  =  first  stress  invariant  =  (0,j  +  (T,,  +  fJ,,) 
b  =  measure  of  distance  between  stress  point  and  surface  (>  1) 

C  =  material  constant  defining  the  projection  center  location 
4  =  bound  size  (i.e.,  intersection  of  bound  with  volumetric  axis) 

J  =bJ 

J  =  second  stress  invariant  s^jS^j 

1  1  1  M 

(2i  —  ^  Cln  —  0  0  ^  _ 

'  2  ^  R  -4X1 

M  =  slope  of  critical  state  line  in  triaxial  space  (assumed  the  same 
in  extension  and  compression) 


R  =  shape  of  ellipsoid. 
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This  function  uses  the  soil  mechanics  sign  convention  of  compression  as  positive  and 
tension  as  negative.  The  surface  is  shown  in  Figure  3.1-1. 

3.2  Elastic  Calculations 

For  stress  paths  starting  within  the  elastic  nucleus,  a  closed  form  solution  can  be 
found  to  predict  the  intersection  of  the  path  with  the  elastic  nucleus.  The  beginning  and  end 
of  the  step  are  denoted  by  4  and  4,^^  respectively.  Therefore,  =  4  is  used  for  the 
beginning  of  the  step  and  for  the  end  of  the  step. 

To  find  the  results  of  an  increment  of  strain  with  a  beginning  value  and  an  assumed 
proportional  strain  history  for  r  =  0  ^  1 ,  a  relative  time  scale  is  used  since  viscous  effects 
are  not  present.  The  strains  can  be  written  as: 

The  rates  of  the  total,  deviatoric  and  volumetric  strains  are: 

The  Bounding  Surface  Plasticity  model  for  clays  uses  the  log-linear  relationship  for 
the  volumetric  stress  versus  the  void  ratio  for  volumetric  stresses  greater  than  the  transition 
stress  (/,)  and  a  linear-linear  relationship  for  stresses  less  than  I,.  This  is  shown  in  Figure 
3.2-1  and  discussed  more  thoroughly  in  Section  1.3. 

Now  at  some  point,  t  in  the  interval  Equation  1 .3-8  can  be  written  in  the  form: 


(/>/,) 

(3.2-3a) 

I  =  /j/,(e' -»:)+/, 

(/</;) 

(3.2-3b) 

It  is  assumed  that  the  path  from  time  4  to  t  will  be  entirely  elastic 
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(i.e..  O'"  -61=  A6‘'t ).  Because  this  is  an  elastic  step  (AO"  =  AS),  however,  the  “e” 
notation  will  be  used  to  emphasize  of  the  nature  of  the  elastic  step.  The  volumetric 
expressions  can  be  rewritten  in  terms  of  the  strain  increment  as: 

I  =  (/>/,)  (3.2-4a) 

1  =  P  l,Ae‘^t  +  I,  (/</,)  (3.2-4b) 

The  tangent  bulk  modulus  is  a  function  of  the  pressure  only  and  can  be  found  by 
taking  the  derivative  of  Equation  3.2-3  with  respect  to  the  volumetric  strain: 

B  =  =  (/>/,)  (3.2-5a) 

3dd^  3^ 

B=  jP  (/</,)  (3.2-5b) 

Writing  the  modulus  in  terms  of  the  behavior  at  the  beginning  of  the  step,  the 
following  is  obtained: 

B  =  (/>/,)  (3.2-6a) 


B=  B,  (/</,)  (3.2-6b) 

where  =  tangent  bulk  modulus  at  the  beginning  of  the  step. 

The  shear  modulus  at  time  t  can  be  expressed  in  terms  of  the  bulk  modulus  as: 

G  =  r]B 


where 


3(1 -2t;) 
2{l  +  v)  ■ 


(3.2-7) 
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Writing  this  in  terms  of  the  value  at  the  beginning  of  the  step,  the  following  is 
obtained: 


G  =  (/>/,)  (3.2-8a) 

G  =  G,  (/</,)  (3.2-8b) 

where  G^  =  shear  modulus  at  the  beginning  of  the  step. 

The  relationship  between  the  deviatoric  stress  rate  and  the  deviatoric  strain  rate  is 
given  as: 


Sy=2Gel  (3.2-9) 

Substituting  Equation  3.2-8  into  3.2-9,  the  deviatoric  stress  rate  can  be  written  as; 
s,i=2a,e^'el  U>I,)  (3.2-lOa) 

i„„,=2G,e-  (/</,)  (3.2- 10b) 

Integrating  the  above  expressions  from  0  to  r  yields: 

s,^{t)  =  2G,Ae;f{t)  +  Sy^  (3.2-11) 

where  f{t)  =  (/>/,) 

f{t)  =  t  (/</,) 

Using  the  definition  for/(0  in  Equation  3.2-1 1,  Equation  3.2-4  can  be  rewritten  as; 
I  =  pA8-f(t)g{l)+J,  (3.2-12) 

where  g{I)  =  h  (/>/,) 


«(/)  =  /, 


(/</,) 
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The  following  definitions  will  be  required  for  the  development  of  this  method: 


=  {pAe'fmijf  +  2pAe- f{t)g{i)  /,  -h  II 

(3.2-13) 

l  =  bI  +  {\-b)I^ 

(3.2-14) 

T-  =  {bif  +  2b{\ -  b)IJ  +  (1  -  bf 

(3.2-15) 

(3.2-16) 

11 

(3.2-17) 

where 

Xi 

Now  the  value  of  t  =  t*  is  defined  where  the  stress  path  intersects  the  elastic  surface. 
The  equation  for  the  elastic  surface  is  given  by  Equation  3.1-1  using  the  value  of  b  that 
defines  the  size  of  the  elastic  surface  in  Equation  1.6-5).  Replacing  b  with  and 
t  with  t*  in  Equations  3.2-13  through  3.2-17  and  substituting  these  equations  into  Equation 
3.1-1  gives: 

F  =  9i/1»')  +  <7./(»')+<!o=0  (3.2-18) 

where  =— (26„„„  Gff  Xt 

=  -i>L^2cai  +  2irL„i/3Ae-«(/)  h 

a, 


qo  =  F  evaluated  at  time 


The  above  equation  can  be  solved  for/  if  =/(t*))  using  the  quadratic  formula: 


(3.2-19) 

2^2 

Possibilities  for  these  roots  are  shown  in  Figure  3.2-2.  For  a  strain  increment  in  a 
positive  direction,  the  positive  root  is  needed  (shown  in  Figure  3.2-2a);  hence  the  +  sign  is 
used  in  Equation  3.2-19.  For  a  strain  increment  in  a  negative  direction,  the  roots  would 
reverse;  hence,  the  +  sign  would  still  be  used.  For  stress  points  close  to  the  elastic  surface 
and  a  strain  increment  in  a  positive  direction,  the  root  is  near  zero  (as  shown  Figure  3.2-2b) 
and  still  requires  using  the  +  sign  in  Equation  3.2-19.  It  is  important  to  note  that  round-off 
error  could  actually  cause  the  root  to  become  complex  (i.e.,  the  term  in  the  radical  may  be 
negative).  In  that  case,  the  root  is  then  taken  to  be  zero.  For  stress  points  close  to  the 
elastic  surface  and  a  strain  increment  in  a  negative  direction,  the  positive  root  is  needed 
(shown  in  Figure  3.2-2c)  thus  requiring  the  +  sign.  When  the  direction  is  tangent  to  the 
surface  (as  shown  in  Figure  3.2-2d)  two  zero  roots  are  the  result.  In  all  cases,  however, 
the  +  sign  is  used  before  the  radical  in  order  to  ensure  a  positive  root. 

At  this  point,  the  length  of  the  elastic  step  (s)  is  sought.  The  value  of  s  will  be 
equal  to  the  value  of  t*  where  the  stress  path  intersects  the  elastic  surface,  or  the  upper  limit 
of  t  =  1  when  the  step  is  entirely  elastic.  Thus,  two  questions  must  be  considered  in 
assigning  the  elastic  step  size: 

1)  Is  the  indicated  step  greater  than  1?  (indicating  that  the  strain  step  is 
entirely  elastic,  as  in  Equation  3.2-1) 

2)  Is  the  transition  stress  crossed? 

One  approach  is  to  evaluate  7  (from  Equation  3.2-12)  and  test  for  the  crossing  of  the 
transition  stress.  However  when  the  indicated  step  size  is  very  large,  numerical  problems 
could  occur,  especially  for  7  <  0.  Therefore,  the  step  size  tests  are  done  on  the  variable/ 
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where/is  defined  in  Equation  3.2-1 1.  The  limit  on  the  value  of /for  ^  <  1  is; 


e^^-\ 

'  pAO 

(/>/) 

(3.2-20a) 

/i=l 

U<Ii) 

(3.2-20b) 

To  determine  whether 

the  stress  path  has  crossed  the  transition  stress  (/),  the  value 

of/ when  the  transition  stress 

line  is  crossed  is  found  as  follows; 

(/>/) 

(3. 2-2  la) 

f- 

‘  p  AO  7, 

(/</) 

(3.2-21b) 

If  the  path  has  crossed  the  transition  stress,  the  elastic  calculation  is  done  as  two 
steps  (e.g.,  the  first  step  is  up  to  /  and  the  second  is  beyond). 

The  quantity/*  is  first  found  using  Equation  3.2-19,  then  compared  against 
Equations  3.2-20  and  3.2-21.  The  smallest  value  of/is  taken: 

/m,„  =  min  [/*  /,  /  ]  (3.2-22) 

The  elastic  step  size  (^)  is  determined  from  Equation  3.2-1 1; 

5  =  -^ln(/3Ae7„,.+l)  (/>/,)  (3.2-23a) 

J  =  (/</,)  (3.2-23b) 

Once  the  value  of  s  is  found,  t  is  set  to  s  and  this  result  is  substituted  into  Equation 
3.2. 12  to  find  the  value  of  the  volumetric  stress  When  7  >  7,  and  Ad  =  0,  the 
equations  simplify  to: 

^elas 


(3.2-24) 


The  stresses  at  the  end  of  the  elastic  step  are  determined  by; 


=<"«.  +2G,A<;;  f(s)+iAIS,  (3.2-25) 

where  AI  =  p  AG‘  f(s)  g{I) . 

Once  the  stresses  are  calculated,  the  loading  surface  measure  (b)  is  updated  along 

with  the  stresses.  If  5  <  1.0  then  the  remainder  of  the  step  (i.e.,  r*  <  t  <  1)  is  plastic.  The 
updated  values  etc.)  are  now  the  beginning  values  for  the  start  of  the  plastic  step  (or 
the  beginning  values  of  the  second  elastic  step  if  the  transition  stress  is  crossed).  The  flow 
chart  for  the  elastic  evaluation  is  shown  in  Figure  3.2-3. 

3.3  Elastic  Contribution  to  the  Global  Jacobian 


Once  the  length  of  the  elastic  step  (5)  is  found,  the  elastic  contribution  to  the  global 
must  be  calculated.  Recalling  the  development  of  the  elastic  step 


Jacobian  matrix 


^  d<7  ^ 

dASr, 

V  "y 


in  Section  3.2,  the  elastic  Jacobian  is  also  evaluated  in  terms  of  deviatoric  and  volumetric 
components.  The  decomposition  of  the  stress  tensor  is  given  as: 


(3.3-1) 


where  <jy  =  total  stress 

Sy  =  deviatoric  stress 
I  =  volumetric  stress  ((J,j  +  (J22  +  CTj^)- 


Differentiating  the  stress  at  s  with  respect  to  the  increment  of  strain  results  in: 
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The  volumetric  portion  can  be  evaluated  by  differentiating  Equation  3.2-4  with 
respect  to  the  strain  increment  and  setting  t  =  s: 


-^  =  pIsS,,  (/>/;)  (3.3-3a) 

dAe^ 


A 

dAe„ 


=  Iil3s5,i 


(/</,)  (3.3-3b) 


The  deviatoric  portion  is  evaluated  by  differentiating  Equation  3.2-1 1  with  respect 
to  the  strain  increment: 


ds^. 


dAe, 


■  =  2G, 


dAel 


dAe,,  l3Ae‘ 


■1  dAd‘ 

— I-- 


dAe„  A6‘ 


se 


PAS^s 


Pass's 


-I 


pAe‘ 


jj 


(3.3-4a) 


ds,j  dAet 

dAe„  dAe,, 


(/</,)  (3.3-4b) 


For  stress  paths  that  cross  through  the  transition  stress  (I,),  the  values  on  both  sides 
oil,  contribute  to  the  global  Jacobian.  In  order  to  develop  their  respective  contributions, 
consider  the  current  stress  which  is  a  function  of  the  previous  stress  and  strain  increment: 


0'y=<Ty(o-*4,Ae^)  (3.3-5) 

Now  if  two  (or  more)  steps  are  taken  (i.e.,  CT*  ->  <7^^  ),  the  final  stress  can  be 

written  as  a  function  of  the  intermediate  stress  and  a  portion  of  the  strain  increment: 


a 


(3.3-6) 


where  s,  =  length  of  first  portion  of  the  step 

S2  =  length  of  second  portion  of  the  step. 
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Differentiating  Equation  3.3-6  with  respect  to  the  strain  increment  and  applying  the 
chain  rule  results  in  the  Jacobian  matrix: 


dAs,,  ‘  dAe, 


d(J 

dASr 


-f 


Ae 


qr 


d{s2-S^) 

dASki 


(3.3-7) 


where  S,- „  =  -- — 


For  simplicity,  the  last  term  in  Equation  3.3-7  (the  term  in  brackets  { })  has  been  neglected. 

Additional  derivatives  that  will  be  required  later  in  the  calculation  of  the  plastic 
contributions  to  the  global  Jacobian  are: 

1)  derivatives  with  respect  to  the  increment  of  strain  of  the  square  of  the  second 
stress  invariant  {f),  and 

2)  the  distance  from  the  current  stress  point  to  the  Bounding  Surface  (b)  with 
respect  to  the  increment  of  strain. 

Differentiating  the  square  of  the  second  stress  invariant  (as  defined  in  Equation  3.2-16) 
results  in: 

=  (3.3-7; 


The  distance  measure  {b)  is  found  by  substituting  Equations  3.2-14  through  -18 
into  the  Bounding  Surface  function  (Equation  3.1-1)  and  solving  the  resulting  quadratic 
equation.  The  equation  for  b  is: 

_ 

a. 


b  =  I. 


(3.3-8) 


78 


For  b  to  be  positive  (and  knowing  that  Uj,  02  and  J  are  always  positive),  the  +  sign 
is  used  before  the  radical.  The  derivative  of  b  with  respect  to  the  increment  in  strain  is 
given,  noting  that  /„  is  fixed  for  an  elastic  step: 


db 

dAs. 


-1 


a 


a, 


(c'-p) 


dl  dJ 


dA£f.i 


a,  dASi.1  j 


a. 


a/' 


dl 


BAs, 


ki ) 


^  dJ  1  BJ^ 

where  — - = - - - 

BAe^,  2J  <7Ae„ 

3.4  Plastic  Calculations 


(3.3-9) 


The  plastic  part  of  the  increment  begins  at  the  end  of  the  elastic  step.  The  stresses 
and  strains  are  updated  after  the  elastic  step.  The  beginning  of  the  plastic  step  is  designated 
with  a  subscript  b  and  treated  as  a  new  step  with: 

t^=0,  <r- Ae-jS,  Ae^  <— (1  — j)Aey  ,  AP<— (1  — 5)A0  (3.4-1) 

To  assure  accuracy,  the  plastic  step  will  allow  substepping.  The  strains  are  sub¬ 
divided  into  N  uniform  substeps  as  follows: 

Ae 

Ae,  (3.4-2) 

‘J.  ^ 

The  notation  used  throughout  this  section  is  1( tj  =  for  the  beginning  of  the  step, 
I(tJ  =  for  the  beginning  of  the  substep,  and  for  the  end  of  the  substep.  The 

plastic  volumetric  strain  can  be  written  in  terms  of  the  total  strain  and  elastic  strain  as: 


(3.4-3) 


79 


Rewriting  Equation  3.4-3  in  terms  of  the  strain  increment  gives; 

e'’it)=e,+Ae{t-t,)-ei-Ae^{t-t,)  (3.4-4) 

Solving  Equation  3.2-4  in  terms  of  the  elastic  strain  increment  and  substituting  into 
Equation  3.4-4  yields; 


d’’{t)=e,+Ae{t-t,)--in 


\  h  J 


-e: 


(/>/,)  (3.4-5a) 


-et 


(!<!,)  (3.4-5b) 


Differentiating  these  equations  with  respect  to  time  (noting  that  the  beginning 
values  are  fixed)  results  in  the  rate  form  relationship  between  the  total,  elastic  and  plastic 
volumetric  strains.  This  relationship  is  written  as  follows; 


1  /, 


n+1 


P  I, 


rt+1 


(/>/,)  (3.4-6a) 


AS-  —  — 

P  h 


(/</;)  (3.4-6b) 


The  associative  flow  mle  for  the  volumetric  plastic  strain  rate  in  the  Bounding 
Surface  formulation  is  given  as; 


QP  =:3y 


n+\ 


(3.4-7) 


where  =  plasticity  parameter  (loading  index). 
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Using  Equation  3.1-1,  the  derivative  of  F  with  respect  to  I ,  taken  at  can  be 


expressed  as: 


^  =2(7 


■pK.) 


(3.4-8) 


where  7..,  =i„i(2..i 

Combining  Equations  3.4-6,  3.4-7  and  3.4-8  yields  the  volumetric  differential 
equation  that  must  be  satisfied  during  plastic  deformation  (where 


67„.,(L,-py„„J  =  Ae--^ 


(/>/,)  (3.4-9a) 


(/</,)  (3.4-9b) 


At  this  point,  the  unknowns  in  Equation  3.4-9  include  the  volumetric  stress,  the 
volumetric  stress  rate,  the  loading  surface  measure,  the  plastic  parameter  and  the  bound  size 

(/„+],  /„+! ,  ^„+i ,  7„+i .  ,, .  respectively).  The  evolution  of  the  Bounding  Surface  size  for 

clays  (4)  is  described  in  Section  1.5  [Kaliakin,  1985;  Dafalias  and  Herrmann,  1986]. 
Equation  1.5-8  written  in  terms  of  the  strain  at  the  beginning  of  the  step  is  given  as: 


L  =Le 


(3.4-10) 


Noting  that  A0„4i  =  ^n+i  “  4  =  Equation  3.4-5  is  now  substituted  into 

Equation  3.4-10  resulting  in: 


^o„^,  ^Oi,  ^ 


P  K 


(/>/,)  (3.4-1  la) 


4..,  =4/ 


1  [ 


n  I,  ) 


(3.4-1  lb) 


For  the  special  case  of  the  volumetric  stress  (/)  crossing  the  transition  stress  (/,),  the 
equations  can  be  modified  to  incorporate  the  transition.  This  case  requires  adding  the 
known  distance  (from  the  current  stress  point  to  the  transition  stress)  to  the  unknown 
distance  (the  term  in  parentheses  in  Equation  3.4-1 1).  For  stress  points  originating  in  the 
log-linear  region  (Equation  3.4-1  la),  the  known  distance  is  in  the  log-linear  portion  and  the 
unknown  distance  is  in  the  linear-linear  portion.  For  stress  points  originating  in  the  linear- 
linear  region,  the  reverse  is  true.  The  modification  is  given  as: 


InlA+iiiil-A 

h) 


(  4  >/,and /„,,</,)  (3.4-12a) 


L,  =  ^ 


( 4  </,and /„,,>/,)  (3.4-12b) 


The  unknown  volumetric  stress  rate  [4^, )  can  be  approximated  with  a  backward 
differences  formula: 


(3.4-13) 


Using  Equations  3.4-1 1  through  -13,  the  unknowns  in  Equation  3.4-9  are  now 
reduced  to  just  ,  4^1  and  . 

Similar  to  the  volumetric  strains,  the  deviatoric  total,  elastic  and  plastic  strain  rates 
at  the  end  of  the  step  are  related  via: 


■  e‘.. 

‘Jn*\ 


(3.4-14) 


The  elastic  deviatoric  stress-strain  rate  relationship  is  given  as: 
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It  is  important  to  note  that  the  shear  modulus  (G„+/  defined  in  Equation  3.2-8)  is  a 
function  of  the  bulk  modulus  which,  in  turn,  is  a  function  of  the  volumetric  stress 
The  associative  flow  rule  for  the  deviatoric  plastic  strain  rates  in  the  Bounding 
Surface  formulation  is  given  as: 


=  y 

ij„^\  *  n 


-^1 


(3.4-16) 


The  derivative  of  the  Bounding  Surface  function  (F)  with  respect  to  the  deviatoric 
stresses  can  be  expressed  in  terms  of  the  invariants  as: 


^  dF^ 


(3.4-17) 


r 


where 


dT  a; 


V  ^  A+i 


An  approximation  for  of  Equation  3.4-15  is  made  by  a  backward  difference 

formula  resulting  in; 


At 


(3.4-18) 


‘n+l  ‘n 


]_ 

N' 


Substituting  Equations  3.4-15  though  -18  into  the  strain  rate  equation  (Equation 
3.4-14),  and  recalling  that  the  deviatoric  strain  rate  is  Ae^  =  e^,  yields: 


Ae..  =  ^ — — 
2G„,.At 


(3.4-19) 
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The  expression  can  be  solved  for  the  deviatoric  stresses  at  the  end  of  the  step  as: 


'2,  Cx  ,  ..  +  S- 

rt+I  y  tJn 


(3.4-20) 


a. 


n+\ 


The  square  of  the  second  invariant  is: 


j'^  ~  ~  s  s 


(3.4-21) 


Substituting  the  deviatoric  stress  defined  in  Equation  3.4-20  into  Equation  3.4-21 
results  in  the  deviatoric  equation  that  must  be  satisfied  during  plastic  deformation: 


'nt\ 


(3.4-22) 


where  Xi=j^e,jAey 


The  unknowns  in  these  equations  include  the  new  shear  modulus  (which  is  a 
function  of  the  volumetric  stress),  second  stress  invariant,  loading  surface  measure,  and  the 
plastic  parameter  ( and  7„^j ,  respectively).  The  square  of  the  second  invariant 
can  be  solved  for  by  rewriting  the  Bounding  Surface  equation  (Equation  3.1-1)  as: 


4.1  = 


1 


^2^n+\ 


(3.4-23) 


This  reduces  the  number  of  unknowns  by  defining  ,  but  introduces  . 
Recall,  however,  that  Equation  3.4-1 1  relates  to  thereby  leaving  the  same  three 
unknowns  as  in  the  volumetric  equation  (Equation  3.4-9):  and  7„^.,. 


At  this  point,  since  there  are  two  equations  (Equation  3.4-9  and  Equation  3.4-22) 
and  three  unknowns,  a  third  equation  is  required.  The  obvious  choice  is  an  equation  that 
relates  the  plasticity  parameter  to  the  amount  of  plasticity  deformation  occurring  within  the 
bound.  The  fundamental  equation  of  the  Bounding  Surface  Plasticity  concept  [Kaliakin, 
1985;  Dafalias  and  Herrmann,  1986]  relates  the  plasticity  to  the  current  stress  state. 
Equation  1.2-4  can  be  rewritten  in  terms  of  the  invariants  as: 


r„.i  = 


(Of' 


4+1  + 


dF'^ 

V  )n+l 


4+i4+i 


(3.4-24) 


The  definition  of  the  hardening  modulus  Jis  given  in  Equation  1.6-4 

[Kaliakin,  1985;  Dafalias  and  Herrmann,  1986].  The  volumetric  stress  rate  (/„+,)  is 

approximated  in  Equation  3.4-13,  and  the  second  invariant  rate  evaluated  by 

taking  the  derivative  of  Equation  3.4-23  with  respect  to  the  current  time 


=  _ 


2  -  -P/,,.,)(L,-p4.,) 


'n-¥\ 


^2^n+l 


(3.4-25) 


The  time  derivatives  of  the  Bounding  Surface  size  are  achieved  by  differentiating 
Equations  3.4-11: 


‘.J 


^  1  7  ^ 


(I>1,)  (3.4-26a) 


n+\  J 


4-+I  4_+i  ^ 


^  1  7  ^ 

A0-  — -2±L 

P  4 


(7<7,)  (3.4-26b) 


For  the  special  case  of  the  volumetric  stress  crossing  the  transition  stress  (Equation 
3.4-12),  the  “a”  and  “b”  equations  are  reversed  (i.e..  Equation  3.4-25a  is  used  for  7  <  7, 
and  Equation  3.4-25b  is  used  for  7  >  I). 
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The  rate  of  the  loading  surface  measure  (4+i )  is  approximated  with  a  backward 
difference  formula  as: 


^  ^  .^2±L_A  (3.4-27) 

At 

Using  Equations  3.4-25  through  -27,  the  number  of  unknowns  for  Equation  3.4-24 
is  reduced  to  the  same  set  of  unknowns  as  in  the  volumetric  and  deviatoric  equations 
(Equations  3.4-9  and  3.4-22),  and  7„^,.  The  equation  is  purposely  written  so  that 

the  hardening  modulus  [k^  _ )  appears  in  the  denominator.  This  is  necessary  because  as 

the  stress  state  approaches  the  elastic  bound,  the  hardening  modulus  approaches  infinity  at 
a  rapid  rate.  When  the  hardening  modulus  is  in  the  denominator,  the  plasticity  parameter 
approaches  zero  as  the  elastic  bound  is  reached. 

Since  the  equations  are  extremely  complex,  the  solution  for  the  three  unknowns  is 
achieved  with  Newton-Raphson  iteration.  The  approach  used  in  this  study  is  to  treat  7  and 
b  as  independent  variables  in  the  volumetric  and  deviatoric  equations  (Equations  3.4-9, 
3.4-22).  The  Bounding  Surface  equation  (Equation  3.4-24)  is  solved  for  the  plasticity 

parameter  (Y).  Another  approach  that  was  used  for  the  Cam-Clay  model  [Herrmann,  1997] 

is  discussed  in  Appendix  B.  The  residuals  are  defined  as  the  error  in  the  governing 
equations  (Equations  3.4-9  and  3.4-22)  and  are  given  as: 

(/>/,)  (3.4-28a) 


(/</,)  (3.4-28b) 
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( 


R,= 


l  +  2G„^,Ar— 

J 


JL  -4G„^.,Ar^Z. -2G„,Atx,-J:  (3.4-29) 


The  Newton-Raphson  method  iterates  on  the  independent  variables  in  order  to  find 
the  roots  of  the  residuals.  The  general  one-dimensional  form  is  given  as: 


■^,+1  = 


/U) 


(3.4-30) 


fW 

where  i  =  iteration  counter 

=  improved  value  of  x 
x-  =  previous  value  of  x 
f(x.)  =  residual  (function  of  x) 
f(x)  =  derivative  of  residual  with  respect  to  x. 

The  two-dimensional  form  involves  matrices  and  matrix  inversion  and  is  given  as 


^'1 

i+l,«+l  i.n+I 


where  'F 


,  1 

[/?,} 

1 

L  2  J 

■(9R, 

dl 

db 

dR^ 

dR^  ■ 

L  dl 

db  ^ 

(3.4-31) 


The  local  Jacobian  (^  is  defined  by  taking  derivatives  of  the  residuals  with  respect 

to  the  independent  variables  I  and  b  (note  that  when  taking  derivatives  with  respect  to  one 
independent  variable,  the  other  independent  variable  is  held  constant).  Dropping  the 


iteration  counter  (i)  for  clarity,  these  derivatives  are  given  as: 


V 


=  6r, 


n+i 


-P\ 


^  dl,  ^ 


V  dl  j 


y  n+\ 


(/>/,)(3.4-32a) 


^1, 


n+\ 


J; _ ^n+\ 

^n+lj 


V  y„+i 


=  6y, 


/i+I 


-P\ 


— 1 
V  A+i  J 


(/</,)(3.4-32b) 


+- 


pAtl, 


n+l 


V  y„+i 


=  6/, 


«+l 


V^^A+i 


dbJn^X 


(3.4-33) 


V  dl  )„+i 
/ 


=  4At^h 


n+l 


a 


l  +  7„„ 


l  +  2G„,,At-^7„,,Z?, 


Vra/M 


V 


V  y„+i 


8G„,iAr 


V  a/  y„^, 


+  G. 


n  +  { 


^  dy^ 

\dl  A+i; 


7, 


n+1 


V  a/  y„^i 


-2Ar;i^, 


7^^ 

a/  j„+i 


(3.4-34) 


=4Ar^2.G 


db 


«+I 


l  +  2G„„Ar^7„^,^„,i  7„.,+^„+i 


f  ay 


a, 


vaay„+,; 


7, 


n+1 


+ 


CL'^ 


l4“2G„+jA/ 

a, 


2^ 


V 


(3.4-35) 


The  derivatives  of  the  plasticity  parameter  are  given  as: 


l  SI  L,  Kl 

J_ 

K„ 


1 


vv  a/  7„+i 


K+i  + 


^  dF^ 

A+I 


Vjr-  \ 


b„+ij„+i 


7V 


a/ 


/n+l 


a/a/ 


2„+i  + 


raF^  fa/^ 


\dl 


a/ 


+ 


,a7' 


^rt+i 


/rt+i 


V  Jn^lJ 


(3.4-36) 
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j_ 

K. 


Y— 1 

Jn  +  l 


4+1  + 


^  dF'^ 
ydJ^  A+i 


4i+i4+i 


A  A+i 


4+1  + 


^  dF^ 

ydT  J„^, 


4^1  + 


Jn+l 


^/i+l 


(3.4-37) 


The  local  Jacobian  introduces  a  large  number  of  derivatives.  These  are  summarized 
below  and  are  grouped  by  terms.  Some  of  the  terms  were  defined  previously  (and  hence, 
their  original  equation  numbers  are  used),  but  they  are  grouped  here  for  clarity.  The 
pertinent  derivatives  of  the  volumetric  stress  (/)  are  given  as: 


r  _  4+1  4 

^n+1  ~  ' 


At 


(3.4-13) 


dl 


«+i 


J_ 

A? 


(3.4-38) 


4+1  ==  4+1  (4+1  -  <^4„., )+ 


(3.4-39) 


dl 

dl 


n+\ 


K., )c^\ 


n+I 


(3.4-40) 


db 


=  4+,-C-4. 


n+l 


(3.4-41) 


Li  =  4+1  (4+1  -  c  ) + 4^]  (4+,  -  c  _ ) + c 


(3.4-42) 


dl 

dl 


=  b  , 

fi-ci 

\ 

h  ,1 

'di 

4 

+  C^ 

^n+] 

n+] 

X 

n-fl  J 

n+l 
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,  * 

n+l 

n+l  J 

dl 

n4-l 


(3.4-43) 


dl 

db 


n4-l 


db 

db 


1-C 


«+l  \ 


dl 


+(4+,-C4j 


n+l . 


(3.4-44) 
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The  derivatives  of  the  other  independent  variable  {b),  the  loading  surface  measure, 
are; 


b 


«+i 


(3.4-27) 


* 

db 


n+1 


_1_ 

Ar 


(3.4-45) 


The  derivation  of  the  shear  modulus  (Equation  3.2-8)  shows  that  for  the  log-linear 
range,  the  modulus  is  a  function  only  of  the  volumetric  stress.  The  derivative  of  the  shear 
modulus  with  respect  to  the  volumetric  stress  is  given  as: 


V  Jn+l  h 


(/>/,)  (3.4-46a) 


^  dG^ 

V  dl  Jn^X 


=  0 


(/</,)  (3.4-46b) 


The  derivatives  of  the  Bounding  Surface  function  (JF)  are: 


V  )  n+\ 


(3.4-47) 


^  dF^ 


V  Jn+l 


(3.4-48) 


=  2 


-p\ 


V  dl  J„^i  J 


(3.4-49) 


^dldb  j 


=  2 


n+l 


y^^Jn^x 


(3.4-50) 
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The  second  invariant  {f)  and  its  derivatives  are: 


r2  _  1  [.2 


+  1  i2  P^o 

^2^n  +  \  ^ 


(3.4-23) 


— 1  =^—  I  f— 1  -^\(Ki-p^  )||  — I  -pf— 1  (3.4-51) 


2/„^. 

2 

V. 

^2^n+\ 

(3.4-52) 


jl  _  _  /7+1  jl  1 _ 

'^n+\  ”  7  ^ 

K+\  ^2^«+l 


^l{^n+l  ){^n+l  (3.4-25) 


'il  /  4/ 

Dr  " 

V  oi 


^  I  2^ 

^«4l  I  J„41 


^ydiLy" 


In.^-Pl,  (3-4-53) 


(J  r  ) 

-«.(44,-p/,j  ^  -P^ 


2_  2  ^2^,2  _^f^l 

A  7  '^n  +  1  T  2  ''«  +  i  ,  -)7 

ArZ7„,, 

“  72  ~  ^l{jn+\~  Ph,,^{jn+\~  PK,^^ 

“2^n+l  '- 


2^1  I 

hbl.i 


T.  (4  +  1  py^,*,  )‘*'('^n+l  P^o,,^ 


(3.4-54) 


The  Bounding  Surface  size  (/„)  derivatives  are: 


h>  -  4  ^ 

^4+1  *4 


f  A0r„„-- 

T  p[  1, 


Z'  \ 
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V  cfi  y„+,  p 


iS4 


(/>/,)  (3.4-1  la) 


(/</,)  (3.4-1  lb) 


(/>/,)  (3.4-55a) 


^91.  ^ 


V  y„+i 


14^ 

P  I, 


(/</,)  (3.4-55b) 


4„+i  4,,.,  *3 


1  7  ^ 

Ae- 


V 


(3/, 


n+I  y 


(/>/;)  (3.4-56a) 


4„*)  4„*i  ^ 


1  7  ^ 

Ad- 


V 


;  y 


(7<7,)  (3.4-56b) 


V  y„., 


4)  , 

V  4+1 


17^ 
Ad 


V 


iS/. 


n+i  y 


-—7 

j8 


1  7  ^ 

— ^ - (7>7,)  (3.4-57a) 

vA^4+i  4+1/ 


^^7  ^ 


V  4+1  V  «/  y„^, 


^7 


Ae--!-  4+' 


P  h 


-i/ 


/  y 


P  At  7; 


(7<7,)  (3.4-57b) 


For  the  special  case  of  the  volumetric  stress  crossing  the  transition  stress.  Equation 
3.4-12  is  used  for  7^  and  the  time  derivatives  are  reversed  (i.e.,  Equation  3.4-56a  is  used 
for  7  <  7, and  3.4-56b  is  used  for  7  >  7,). 

The  Newton-Raphson  method  is  set  up  using  Equations  3.4-28  and  -29  as  the 
residuals  and  and  4+/  as  the  independent  variables.  The  method  is  begun  by  assigning 
7„^;  and  their  current  values  and  calculating  a  new  estimate  with  Equation  3.4-31.  The 
iteration  continues  until  a  convergence  criterion  is  met.  The  criterion  chosen  for  this  study 
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is  that  the  change  in  residuals  be  smaller  than  a  user-defined  percentage  of  the  absolute 
value  of  the  total  strain  increment.  This  is  given  as: 


1^1 1  + 


(3.4-58) 


where  ||i?||  =  norm  on  the  residuals. 

This  requires  that  the  residuals  are  in  terms  of  strains,  or  at  least  their  magnitude  is 
that  of  strain.  Another  advantage  of  having  the  residuals  near  the  same  magnitude  is  to 
prevent  any  numerical  problems  in  the  matrix  inversion.  Residual  1  (Equation  3.4-28)  is 
already  in  terms  of  strain.  Residual  2,  however,  (Equation  3.4-29)  is  in  terms  of  stress. 
This  residual  is  given  the  magnitude  of  strain  by  dividing  it  and  its  derivatives  by  the  initial 
shear  modulus.  The  new  residual  and  derivatives  are  given  as: 

n*  _  ^2  ^^2  _  ^  ^  ^^2  (3.4-59) 

^  GI  ’  dl  GI  dl  ’  db  GI  db 

In  order  to  add  robustness,  the  plastic  increment  is  uniformly  substepped  as 
described  in  Equation  3.4-2.  First  a  single  step  is  tried,  then  the  strain  increment  is  divided 
by  two  and  the  stress  recalculated.  This  process  continues  until  some  convergence  criterion 
is  met  or  the  maximum  number  of  allowable  substeps  is  exceeded.  In  order  to  compare  the 
various  methods,  a  criterion  based  upon  terms  that  are  present  in  each  of  the  methods  was 
chosen.  The  criterion  used  for  the  substepping  evaluation  is  that  the  change  of  the 
■  Bounding  Surface  size  (/„)  should  be  less  than  a  user-defined  tolerance.  This  is  written  as 
follows: 


THEN:  EXIT  (3.4-60) 


where  I  =  stress  at  end  of  step  and  substep  level,  m 

I  =  stress  at  end  of  step  and  substep  level,  m-1 


93 


m  =  substep  level  (i.e.,  number  of  substeps  =  1,  2,  4,  8,  16,  ...) 

It  =  iteration  counter 

=  maximum  number  of  iterations. 

Another  criterion  considered  involves  the  stress  distance  measure  The  value 

of  is  used  during  the  iteration  as  an  independent  variable.  It  can  also  be  calculated 
independently  knowing  the  stress  state  and  the  Bounding  Surface  size  (Equation  3.3-8).  A 
convergence  criterion  can  be  defined  as  follows: 

{[(*.., or  [/»/<„,]}  THEN:  EXIT  (3,4-61) 

where  =  b  calculated  from  and  ^ . 

This  criterion  is  important  for  stress  points  near  1  =  0,  J=0  and  C  =  0  because  the 
iteration  residuals,  and  thus  the  iteration  criterion  (Equation  3.4-58),  are  near  zero,  but  the 
variable  b  can  vary  significantly. 

An  important  difference  exists  in  the  use  of  the  subscripts  b  and  n.  Most  of  the 
equations  for  the  method  are  defined  in  terms  of  the  beginning  of  the  plastic  increment 
(e.g.,  4).  Terms  defined  at  the  beginning  are  fixed,  thus  simplifying  derivations  (i.e., 
terms  defined  at  time  do  not  contribute  to  derivatives  taken  at  time  The  rate 

equations  for  the  independent  variables  /  and  b  (Equations  3.4-13  and  3.4-27),  however, 
require  information  at  the  beginning  of  each  substep  (i.e.,  4  and  bj.  Therefore,  the 
independent  variables  at  the  end  of  the  substep  must  be  reserved  for  calculating  these  rate 
equations  and  for  eventual  use  in  computing  the  algorithmic  consistent  moduli  (Section 
3.5).  The  algorithm  is  given  in  Box  3.4-1. 


Box  3.4-1.  Reduced  Newton  Algorithm  for  Plastic  Steps. 


1. 

2. 

3. 

4. 

5. 

6. 

7. 

8. 

9. 

10. 

11. 


Solve  for  Increment  in  Variables: 

Calculate  and  Test  Norm;  IF  (||R||  <  Toler,) 

TRUE;  Go  to  step  1 1 

FALSE:  Increment  variables,  Go  to  step  7 

Calculate  Stresses  and  Check  Substeps: 

IF  (k<N)  TRUE:  Go  to  step  4  FALSE:  Go  to  step  12 

Compare  stresses  with  previous  substep  level: 


TRUE:  EXIT 


FALSE:  N  =  2N,  Go  to  step  2 


3.5  Plastic  Contribution  to  the  Global  Jacobian 


Once  the  plastic  step  has  converged,  the  plastic  contribution  to  the  global  Jacobian 


matnx 


must  be  calculated.  The  local  Jacobian  will  be  evaluated  in  terms  of  the 


deviatoric  and  volumetric  components  similar  to  the  elastic  Jacobian.  The  volumetric 
component  of  the  Jacobian  can  be  obtained  by  noting  that  the  residuals  are  a  function  of  the 
independent  variables  I  and  b  which,  in  turn,  are  a  function  of  the  strain  increment.  The 
derivatives  can  be  obtained  by  an  application  of  the  chain  rule:  holding  the  independent 
variables  constant  and  then  adding  the  derivatives  with  respect  to  the  independent  variables. 
At  convergence,  the  residual  is  equal  to  zero  and  the  derivatives  can  be  written  as: 


r 


dR„ 


V  J 


dR 


n+I 


dR,. 


+  ■ 


Expressed  in  matrix  form: 


dl 


^  dl  ^ 


dR,. 


+  ■ 


db 


db 


=  0  (3.5-1) 


dRi  ' 

dl  1 

dA£(.i 

dR, 

>  =:  < 

dAShi 

dR, 

dAs^ 

db 

,  dAS/j 

n+1 

dASki  \ 

^n  +  l  ’^n  +  1 

fOl 

lol 


n+\ 


(3.5-2) 


where  'F 


dR^ 

dR^ 

'W 

db 

dR, 

dl 

db  . 

Note  that  the  matrix  ('P)  has  been  computed  in  the  stress  point  algorithm  (Equation 
3.4-31).  Equation  3.5-2  can  be  rearranged  to  solve  for  the  derivatives  of  the  independent 
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variables  with  respect  to  the  strain  increments.  This  is  given  as: 


(9/  1 

<99?,  ' 

<9Aey 

db 

1 

1! 

dAe^i 

<9/?2 

n+l 

(3.5-3) 


The  derivative  of  the  first  residual  (Equation  3.4-28)  with  respect  to  the  strain 
increment  (holding  the  independent  variables  at  the  end  of  the  step  b^^,)  constant)  is 
given  as: 


SR, 


(/>/,)(3.5-4a) 


6(Li-P4,J: 


-4- 


^  +  1  +  1 


_J _ 

SACu  ,  , 


=  5y 

dAe., 


—  D - 

<9Ae,„ 


1  <9/ 

_ n 

PAtli  dAe^  , 


(/</,)(3.5-4b) 


The  derivative  of  the  second  residual  (Equation  3.4-29)  is  given  as: 


=  4G„,,Ar^^„„[  1  +  2G„,,aA7„,A>i 
+fl  +  2G„„AA7„^i^„,il 

a,  J  <9Ae,,  ,  , 


(3.5-5) 


I..-.- 


•2G„^,Ar-:^ 

dAe^  dAe^ 
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The  derivative  of  the  plasticity  parameter  with  respect  to  the  strain  increment  is 
given  by: 


dj, 


n+\ 


dAe, 


^n+ ! »^n+l 


1  fdF  i  ,  <9F  ^  -2 

dJ  j 


1 


,  dJl 


dJ  dAe. 


K 


Prt  +  \ 

2 


1  dF  dF 


At  dl  dAe^, 


V 


,  ,  3P  -“aAe^ 

^n+l ’^n+1 


4  +  1  ’4i  +  l  / 


F, 


n  +  1 


-P- 


^A£„  ^Ae, 


P„., 


(3.5-6) 


Additional  derivatives  required  for  Equations  3.5-4  through  -6  are  given  as: 


dl 


n+I 


dAe, 


4+1  ’4i+i 
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4+1  ’4i+i 
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SAe,, 
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(3.5-8) 
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dAe, 


(3.5-9) 
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dl 


^  P  -^n+l 


^n+1  ’^fl  +  1  y 


(/>/,)  (3.5-lOa) 
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strain 


<9/„ 


L.,.h„ 


A6»- 1-^211 
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(3.5-11) 
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dAe, 
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2b  bJl, 

K.xK 
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""*■  dAe, 


+  l 


dl 


dAe, 
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-«i(Li-p4 
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1 

dASki 

if 

4  +  1  ’^n  +  !  / 

(3.5-12) 


dx^ 


dAe, 


dx, 


=  e. 


L.,.b.. 


dAe, 


(3.5-13) 


dXi 


dAe, 


L.xA 


-^  =  s  +e 

dAe,, 


^  ds,  ^ 


"*\dAe, 


(3.5-14) 


The  derivative  of  the  deviatoric  stress  at  time  with  respect  to  the  increment  in 


V  dAe, 


in  Equation  3.5-14  is  known  from  the  previous  substep.  If  the  previous 


substep  was  elastic,  the  derivative  is  calculated  via  Equation  3.3-4. 


Note  that  this  derivative  requires  the  derivatives 


and 


calculated  in  Equation  3.5-3.  Also,  a  number  of  the  derivatives  are  taken  at  time  (i.e., 
values  found  at  the  end  of  the  previous  substep  and  saved).  When  substepping  is  used  in 
the  stress  point  iteration,  it  also  must  be  used  in  the  evaluation  of  the  local  Jacobian.  The 
same  substep  level  is  used  for  both  the  calculation  of  stress  and  the  local  Jacobian.  As 
described  in  Section  3.4,  the  independent  variables  1  and  b  are  saved  at  the  end  of  each 
substep.  The  local  Jacobian  is  constructed  by  recalling  these  values  to  calculate  the  current 
substep  derivatives.  These  current  values  are  reserved  and  used  in  the  calculation  of 
the  derivatives  at  time  etc.  For  steps  that  transition  through  the  elastic  zone,  the  first 
substep  requires  values  calculated  at  the  boundary  of  the  elastic  nucleus  (Equations  3.3-4 
and  3.3-9).  At  the  end  of  the  substepping,  both  the  volumetric  and  deviatoric  components 
are  combined  to  provide  the  local  contribution  to  the  global  Jacobian. 
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3.6  Unloading 

For  each  global  step,  it  must  be  determined  whether  the  step  is  loading  (elastic  or 
plastic)  or  unloading  (elastic).  One  common  approach,  especially  useful  when  a  predictor 
stress  is  calculated,  is  to  cross  the  stress  increment  with  the  normal  to  the  function  surface 
(i.e.,  the  loading  index).  A  positive  value  indicates  loading  and  a  negative  value  indicates 
unloading.  The  loading  index  is  given  as; 

L  =  -^A(7  (3.6-1) 

da 

where  L  =  loading  index  (L  >  0,  loading;  L  <  0  unloading). 

Another  approach,  similar  to  the  elastic  case,  is  to  determine  whether  there  are  any 
positive  roots  when  calculating  the  intercept  to  the  loading  surface.  This  approach 
additionally  provides  the  terms  necessary  to  calculate  the  length  of  the  elastic  unloading 
step.  Equation  3.2-18  can  be  rewritten  as: 

^  =  ^2/^(0  +  ^i/(0  +  ^o  =0  (3.6-2) 

where 

c,=^(2G,fz,+(M0g(/)f 

^1 

+b  c-^ 

c,=^2Ga2  +2PAeg(I)I,-2I,0Aeg(I) 

c,=2pAdg{I)(l,-plJ 

qo^b^c^+bc^+Cf, 
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<:s=^(K-Ph...)i.-2K{h-PK..,) 


c^  =  I^-2pIJ^  + 


1 


^ 


a 


\) 


Similar  to  the  elastic  case,  assume  that  t  in  Equation  3.6-2  can  take  on  any  value 
(t**).  The  loading  surface  intercept  is  defined  by  using  the  current  value  ofb  =  b^  and 
solving  for  the  value  of/in  Equation  3.6-2.  Defining/''*  =  the  solution  is: 

2^2 


One  root  is  always  near  zero  (since  the  current  stress  point  is  on  the  loading 
surface)  and  the  other  will  be  either  positive  (indicating  an  intercept  with  the  loading  surface 
and  therefore  unloading)  or  negative  (which  indicates  loading).  Therefore,  the  root  with 
the  larger  absolute  magnitude  is  tested  for  its  sign: 

n  3  /„  =  max  (|/,**|,|/2**|)  (3.6-4) 

where  f**  >  0  =?>  unloading 
f**  <  0  =»  loading. 


If  the  test  determines  that  unloading  is  occurring,  the  length  of  the  unloading  step 
must  be  calculated.  Consider  an  unloading  stress  path  as  shown  in  Figure  3.6-1.  Initially 
the  loading  surface  shrinks  or  unloads  until  the  stress  path  becomes  tangent  with  the 
surface.  At  this  point  the  surface  begins  to  grow  and  loading  begins.  Now  if  b  is  treated 
as  an  unknown  in  Equation  3.6-3,  we  can  determine  at  what  value  of  b  the  two  roots  will 
be  the  same,  thus  defining  the  loading  surface  just  tangent  to  the  stress  path.  The  roots  are 


the  same  when  the  radical  term  in  Equation  3.6-3  goes  to  zero  (i.e.,  -4^29o  ~  ^)-  "This 

yields  an  equation  for  b  using  the  terms  calculated  in  Equation  3.6-2: 


b'\yq^+bq^-\-q^\  =  0  (3.6-5) 

where  q^  =  cl  —  4c,c^ 

^4  =2^203 -4c, C5 

^5  =C2^-4c,C4. 

Noting  that  the  definition  of  b  requires  that  it  be  greater  than  one,  the  two  zero  roots 
are  discarded.  The  quadratic  in  the  square  brackets  is  solved  for  b^^^  where  the  loading 
surface  is  a  minimum; 

min  ^  ^ 

2^5 

Because  the  loading  surface  is  described  in  f  space  (i.e.,  the  square  of  the 
deviatoric  stresses),  one  root  will  be  positive  and  the  other  will  be  negative.  Again, 
because  the  definition  of  b  requires  that  it  be  greater  than  one,  the  positive  root  is  taken. 

The  value  of  b^^^  can  be  substituted  for  b  into  Equation  3.6-2  and  a  new/**  is  determined 
(Equation  3.6-3)  which  yields  the  length  of  the  unloading  step  (5„„,^^)-  The  calculation  of 
the  behavior  of  the  unloading  step  is  the  same  as  described  for  the  elastic  step  calculation 
(Equations  3.2-26).  Solving  for  b  first  is  advantageous  because  if  the  minimum  falls 
within  the  elastic  surface  (i.e.,  b„^  >  b^ias,ic)<  the  elastic  intercept  can  be  easily  calculated  by 
substituting  b  =  The  flow  chart  for  determining  the  elastic  unloading  step  is  shown 

in  Figure  3.6-2. 


Void  Ratio, 


/ 


Figure  3.1-1.  Bounding  Surface  in  Two-Invariant  Space. 


First  Effective  Stress  Invariant,  7 
Figure  3.2-1.  Volumetric  Stress  versus  Void  Ratio  Relationship. 
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(c )  (d ) 


Figure  3.2-2.  Possible  Roots  for  the  Intersections  of  the  Elastic  Surface. 
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Figure  3.2-3.  Elastic  Step  Flow  Chart 


Figure  3.6-1.  Typical  Unloading  Path. 


4.  Example  Calculations 


4.1  Problem  Descriptions 

Comparisons  of  predictions  for  a  number  of  sample  problems  were  made  using  the 
three  methods;  1)  trapezoidal  (Trap),  2)  Closest  Point  (CP),  and  3)  Reduced  Newton 
(RN),  in  order  to  determine  their  behavior.  This  comparison  section  is  divided  into  three 
sections;  1)  comparison  of  the  methods  without  substepping  starting  on  the  bound,  2) 
comparison  of  the  methods  without  substepping  starting  inside  the  bound  and  3) 
comparison  of  the  methods  with  substepping.  Substepping  involves  uniformly 
subdividing  the  strain  increment  within  the  material  model  subroutine  and  solving  for  the 
stress  as  a  cumulative  sum  of  the  incremental  stresses  from  each  substep. 

Section  4.2  compares  the  methods  without  substepping  and  provides  a  general 
description  of  how  each  method  behaves.  Because  the  convergence  criterion  for  the 
iterations  of  each  of  the  methods  is  different,  comparison  of  computational  times  is 
meaningless.  The  methods  are  compared  at  1,  2, 4  and  100  global  steps.  There  is  a 
difference  between  taking  N  substeps  across  the  interval  as  compared  to  taking  N  global 
steps.  The  histories  of  the  various  strain  components  are  generally  not  proportional.  When 

N  global  steps  are  used,  this  nonproportionality  is  taken  into  account  (i.e.,  4  /(O)- 

In  contrast,  the  use  of  N  substeps  approximates  the  strain  histories  as  proportional 
(i.e.,  e^j/£,^i=  constant). 

For  points  that  start  on  the  bound,  the  “exact”  solutions  are  established  by  assigning 
an  arbitrary  change  in  state  and  calculating  (using  numerical  integration  with  a  very  tight 
tolerance)  the  corresponding  strains.  The  strains  for  the  intermediate  steps  are  also 
determined  by  integration.  The  strain  calculations  are  given  in  Appendix  C. 


Section  4.3  considers  stress  points  that  start  within  the  bound.  Since  no  “exact” 
strains  were  calculated,  an  arbitrary  strain  increment  was  chosen.  The  trapezoidal  solution 
for  100  uniform  global  steps  is  assumed  to  be  nearly  exact  for  comparison  purposes. 

Section  4.4  compares  the  methods  with  uniform  substepping  within  the  algorithm. 
A  common  convergence  criterion  is  established  for  the  substepping,  therefore  comparison 
of  computational  times  is  made  possible.  The  criterion  used  is  defined  in  Equation  3.4-59. 
A  set  of  strain  increments  is  chosen  for  each  problem  and  the  algorithm  is  allowed  to 
uniformly  substep  as  required  to  meet  the  substepping  convergence  criterion. 

4.2  Comparison  of  Methods  for  Stress  Points  on  the  Bound  (no 
substepping) 

Several  tests  were  performed  to  compare  of  the  behavior  of  the  methods  for  stress 
points  initially  on  the  bound.  The  Bounding  Surface  parameters  used  for  these  tests  are 
given  in  Table  4.2-1.  Discussion  of  the  parameters  is  given  in  the  references  [Herrmann 
and  Mish,  1983b  and  Kaliakin,  1985]. 


Table  4.2-1.  Bounding  Surface  Parameters  for  Stress  Points  on  the  Bound. 


A 

K 

M 

V 

h 

R 

0.14 

0.05 

1.05 

0.2 

14.7  (psi) 

2.6 

4.2.1  Volumetric  Compression  Test  1 

A  simple  compression  test  along  the  volumetric  axis  was  analyzed  using  the  three 
methods.  The  problem  data  is  given  in  Table  4.2. 1-1.  The  “exact”,  Trap,  CP  and  RN 
results  are  given  on  Tables  4.2.1-2,  4.2.1-3, 4.2.1-4  and  4.2.1-5,  respectively.  The  Trap 
method  provides  reasonable  results  at  one  step  with  a  slight  improvement  as  the  number  of 


steps  increases.  Because  the  CP  and  RN  methods  treat  the  volumetric  behavior  separately, 
they  achieve  the  exact  answer  within  a  few  iterations. 

4.2.2  Volumetric  Tension  Test 

Similar  to  the  compression  test,  a  simple  tension  test  along  the  volumetric  axis  was 
also  performed.  The  problem  data  and  the  “exact”  results  are  given  in  Table  4.2.2- 1  and 

4.2. 2- 2,  respectively.  The  first  attempt  using  the  Trap  method  is  shown  in  Table  4.2.2-3. 
The  solution  is  wrong  and  does  not  improve,  even  with  100  steps.  The  problem  lies  in  the 
implementation  of  the  hardening  rule  (Section  1.5).  In  the  previous  work  [Kaliakin, 

1985],  the  hardening  rule  was  broken  into  a  linear-linear  and  log-linear  section  similar  to 
the  volumetric  stress  behavior  (Section  1.3).  As  strain  developed  for  a  given  stress  path 
(Appendix  C,  Section  C.l),  the  amount  of  tensile  volumetric  strain  that  this  form  of  the 
hardening  rule  could  tolerate  was  limited.  Simplifying  the  hardening  rule  to  the  log-linear 
form  improved  the  performance,  although  a  large  number  of  steps  were  required  to  achieve 
a  reasonable  answer  (see  Table  4.2.2-4).  In  contrast,  the  CP  and  RN  methods  exhibited 
accurate  results  with  one  step  (see  Tables  4. 2. 2-5  and  4.2. 2-6). 

4.2.3  Shear  Test  1 

A  shear  test  was  designed  so  that  the  stress  is  initially  on  the  volumetric  stress  axis 
and  travels  in  the  direction  of  the  second  stress  invariant  (I)  while  keeping  the  first  invariant 
(7)  constant.  The  problem  data  is  given  in  Table  4.2.3-1.  The  results  for  the  Trap,  CP  and 
RN  methods  are  shown  in  terms  of  J  and  the  hardening  (i.e.,  the  bound  size,  IJ  in  Figures 

4.2.3-  1,  4.2. 3-2  and  4.2.3-3,  respectively.  The  Trap  method  shows  the  largest  error  at  the 
larger  strain  increments,  which  can  be  expected  since  the  integration  does  not  explicitly 
satisfy  the  consistency  condition  at  the  end  of  the  step  (i.e.,  maintaining  F  =  0).  The  error 
manifests  itself  in  larger  shear  stresses  and  less  hardening  than  the  “exact”  stress  path. 

This  behavior  indicates  that  the  final  stress  point  is  outside  of  the  Bounding  Surface.  The 


CP  and  RN  methods  show  similar  behavior  because  they  both  are  Newton-Raphson 
methods  that  explicitly  satisfy  the  consistency  condition  at  the  end  of  the  step.  Although 
there  is  some  error  for  the  larger  strain  increments,  the  stress  points  fall  close  to  the  “exact” 
stress  path.  The  convergence  from  below  the  “exact”  path  indicates  points  falling  within 
the  bound. 

4.2.4  Shear  Test  2 

Shear  Test  1  is  now  modified  so  that  the  stress  begins  some  distance  from  the 
volumetric  axis  (i.e.,  J  =  0),  but  still  on  the  bound.  This  data  is  given  in  Table  4.2.4- 1. 
The  Trap  (Figure  4.2.4- 1)  again  shows  larger  shear  stresses  and  less  hardening  than  the 
“exact”  stress  path  for  large  strain  increments.  For  smaller  strain  steps,  however,  the  Trap 
method  exhibits  somewhat  better  accuracy.  The  CP  and  RN  methods  (Figures  4.2.4-2  and 
4.2.4-3)  show  virtually  the  same  behavior  and  converge  from  the  opposite  side  than  the 
Trap  method. 

4.3  Comparison  of  Methods  for  Stress  Points  inside  the  Bound  (no 
substepping) 

Several  tests  were  conducted  for  comparison  of  the  behaviors  of  the  methods  for 
stress  points  initially  inside  the  bound.  The  Bounding  Surface  parameters  used  for  these 
tests  include  those  for  points  on  the  bound  (given  in  Table  4.2-1)  and  the  parameters 
defining  behavior  within  the  bound  (given  in  Table  4.3-1).  Discussion  of  these  parameters 
is  given  in  the  references  [Herrmann  and  Mish,  1983b  and  Kaliakin,  1985]. 


Table  4.3-1 .  Bounding  Surface  Parameters  for  Stress  Points  within  the  Bound. 


P  atm 

C 

m 

a 

w 

14.7  psi 

0.28 

0 

0.02 

4.0 

4.0 

1.2 

5.0 

4.3.1  Volumetric  Compression  Test  2 


A  simple  strain-controlled  compression  test  along  the  volumetric  axis  is  described  in 
Table  4.3. 1-1 .  No  “exact”  solution  was  calculated  (i.e.,  a  stress  path  is  assumed  and  the 
corresponding  strains  calculated)  so  the  Trap  solution  using  100  steps  was  assumed  to  be 
the  actual  answer.  The  Trap,  CP  and  RN  results  are  given  in  Tables  4.3. 1-2, 4.3. 1-3  and 

4.3. 1- 4,  respectively.  All  of  the  methods  provided  reasonable  results  at  one  step  and  a 
slight  improvement  as  the  number  of  steps  increased. 

4.3.2  Shear  Test  3 

A  shear  test  was  analyzed  that  starts  the  stress  on  the  volumetric  axis  and  moves 
mostly  in  the  second  stress  invariant  (J)  direction.  The  problem  data  is  given  in  Table 

4.3.2-  1.  The  results  for  the  Trap,  CP  and  RN  methods  are  shown  in  terms  of  J  and  in 
Figures  4.3.2-1,  4.3. 2-2  and  4.3.2-3,  respectively.  The  Trap  method  shows  greater  shear 
stress  and  less  hardening  for  the  larger  strain  increments.  The  CP  and  RN  methods,  while 
having  significant  error  for  the  large  strain  increments,  show  end  points  that  fall  on  the 
actual  stress  path. 

4.3.3  Shear  Test  4 

Shear  Test  3  is  now  modified  so  that  the  stress  begins  some  distance  from  the 
volumetric  axis  (i.e.,  J  =  0).  This  data  is  given  in  Table  4.3.3- 1.  The  Trap  method  (Figure 

4.3.3-  1)  for  this  test  actually  shows  the  best  behavior  converging  on  the  “exact”  stress  path 
within  two  steps.  The  CP  and  RN  methods  (Figures  4. 2.4-2  and  4.2.4-3)  show  poorer 
accuracy  and  converge  to  the  “exact”  answer  from  the  opposite  side  than  the  Trap  method. 


4.4  Comparison  of  Methods  with  Substepping 


The  Bounding  Surface  parameters  used  for  the  substepping  tests  are  given  in  Tables 
4.2-1  and  4.4-1.  The  only  difference  from  the  solutions  inside  the  bound  (Section  4.3)  is 
that  a  small  elastic  nucleus  has  been  introduced  (e.g.,  =  3.0). 


Table  4.4-1.  Bounding  Surface  Parameters  used  for  Substepping  Tests. 


P  atm 

C 

m 

tic 

a 

w 

14.7  psi 

0.28 

3.0 

0.02 

4.0 

4.0 

1.2 

5.0 

The  substepping  used  for  this  study  is  a  linear  subdivision  of  the  strain  increment 
(i.e.,  dividing  the  strain  increments  by  two)  up  to  a  level  of  32  (i.e.,  2^)  equal  substeps. 
The  substepping  levels  and  number  of  steps  are  given  in  Table  4.4-2. 


Table  4.4-2.  Substep  Levels  versus  Number  of  Substeps. 


Substep  Level,  m 

1 

2 

3 

4 

5 

6 

Number  of  Equal  Substeps 

1 

2 

4 

8 

16 

32 

The  convergence  criterion  used  to  stop  substepping  in  this  study  is  the  percentage 
of  change  in  the  predicted  final  Bounding  Surface  size.  This  was  defined  in  Equation  3.4- 
60  and  is  given  by: 


OR  THEN:  EXIT 


(3.4-60) 


The  methods  iterate  on  the  solution  for  each  strain  level  produced  by  the 
substepping  until  the  iteration  criterion  is  satisfied.  At  the  end  of  each  substepping  process, 
the  Bounding  Surface  size  (IJ  is  compared  to  the  value  found  with  the  previous 


substepping.  When  the  change  in  the  bound  size  is  below  the  specified  tolerance,  the 
results  are  printed.  The  tolerance  used  for  this  study  was  0.001  psi. 

4.4.1  Nearly  Elastic  Compression  and  Shear 

A  test  was  designed  to  start  within  the  elastic  nucleus  and  project  out  into  the  plastic 
realm.  The  initial  stress  is  located  on  the  volumetric  axis  (i.e.,  7=0)  near  the  projection 
center  so  a  large  part  of  the  path  is  elastic.  This  data  is  given  in  Table  4.4. 1-1.  The  stress 
path  in  /  -  7  space  is  shown  in  Figure  4.4. 1-1.  The  “exact”  solution  uses  the  Trap  method 
and  100  linearly  interpolated  global  strain  increments. 

The  Trap  and  RN  methods  show  nearly  identical  results  and  nearly  give  the  “exact” 
solution.  The  CP  method  shows  significant  error  in  the  stress  path.  This  error  is  caused 
by  the  inaccuracy  of  the  rate  equation  for  the  stress  distance  measure  {b)  near  the  elastic 
surface. 

The  computational  time  of  the  methods  were  compared  with  the  codes  compiled 
with  Language  Systems  FORTRAN  on  a  Macintosh  Centris  650.  The  times  are  compared 
on  a  relative  scale  with  the  Trap  method  time  used  as  the  base  (i.e.,  the  Trap  method  is 
always  one).  The  relative  speeds  and  substepping  are  shown  in  Figure  4.4. 1-2.  Although 
the  Trap  and  RN  methods  show  nearly  identical  results,  the  Trap  method  shows  the  best 
behavior  converging  on  the  “exact”  solution  in  less  time  and  fewer  substeps.  This  is  to  be 
expected  since  it  is  a  second  order  method  and  is  most  appropriate  for  nearly  elastic 
calculations.  Also,  the  RN  and  CP  methods  converge  in  the  same  number  of  substeps,  but 
the  CP  method  takes  significantly  longer.  This  is  a  result  of  the  second  order  derivatives 
that  need  to  be  calculated  every  iteration.  In  addition,  each  iteration  of  the  CP  method  must 
solve  nine  simultaneous  equations  while  the  RN  method  only  solves  two.  Another 
difficulty  with  the  CP  method  is  the  use  of  the  rate  equation  for  b  (Equation  1.7-1 1)  which 
is  a  function  of  the  hardening  modulus  {K^).  As  the  hardening  modulus  approaches  the 
elastic  surface  it  rapidly  goes  to  infinity  (see  Equation  1.6-4).  The  behavior  of  the 


hardening  modulus  and  the  b  rate  equation  is  shown  in  Figure  4.4. 1-3.  For  stress  points 
near  the  elastic  surface,  the  use  of  Equation  1.7-1 1  actually  slows  convergence  because  of 
the  large  gradient. 

4.4.2  Plastic  Compression  and  Shear 

A  plasticity  test  was  designed  to  start  near  the  bound  and  project  out  into  the  plastic 
realm.  The  initial  stress  is  located  off  the  volumetric  axis  (i.e.,  J  =  -3.6  psi)  near  the  bound 
so  a  significant  portion  of  the  path  is  on  the  Bounding  Surface.  This  data  is  given  in  Table 
4.4.2- 1.  The  stress  path  in  7  -  7  space  is  shown  in  Figure  4.4.2- 1. 

All  of  the  methods  show  results  close  to  the  “exact”  solution.  The  CP  method 
shows  slightly  more  error  than  the  other  methods. 

The  relative  times  are  shown  in  Figure  4.4.2-2,  and  for  the  Trap  and  RN  methods 
they  are  almost  the  reverse  of  what  was  observed  in  Section  4.4.1.  The  RN  method  shows 
about  a  50%  decrease  from  the  Trap  method.  The  CP  and  RN  methods  were  developed  to 
account  for  plasticity  and  show  improved  performance  by  converging  in  fewer  substeps. 
Again,  because  of  the  calculation  of  the  second  order  derivatives  in  the  CP  method,  the 
relative  time  is  significantly  higher. 

4.4.3  Shear  Softening 

A  shear  test  in  the  softening  realm  was  designed  to  start  on  the  bound.  Volumetric 
strains  were  prescribed  to  keep  the  volumetric  stress  nearly  constant.  This  data  is  given  in 
Table  4.4.3- 1.  The  stress  path  in  7  -  7  space  is  shown  in  Figure  4.4.3- 1. 

All  of  the  methods  show  results  close  to  the  “exact”  solution.  The  CP  and  RN 
methods  show  slightly  more  error  than  the  Trap  method,  but  this  is  small  (note  the  scale  of 
the  7-axis). 

The  computational  times  are  shown  in  Figure  4.4.3-2.  The  CP  and  RN  show  a 
larger  number  of  substeps  for  convergence,  yet  the  RN  method  took  less  time  than  the  Trap 


method  with  fewer  substeps.  Again,  because  of  the  calculation  of  the  second  order 
derivatives  in  the  CP  method,  the  relative  time  is  significantly  higher. 


Table  4. 2.1-1.  Compression  Test  l-Parameters. 


/„  =  8 1  (psi),  e„  -  0.94 

0ij  —  O22  “  ^33  -27  psi,  0]2  0|3  “  O23  —  0 

Ae,,  =  A£22  =  A£j3  =  -2.80025x10  A£,2  =  A£,3  =  A£23  =  0 


J  K 


Table  4.2. 1-2.  Compression  “Exact”  Results. 


1 

h 

b 

91 

91 

1 

Table  4.2. 1-3.  Compression  Test  1 -Trapezoidal  Results. 


Number  of  Steps 

I 

h 

b 

1 

90.9799 

91.0101 

0.9996 

2 

90.9894 

91.0012 

0.9997 

4 

90.9968 

90.9980 

1.0000 

Table  4.2. 1-4.  Compression  Test  1-Closest  Point  Results. 


Number  of  Steps 

/ 

b 

1 

91.0000 

91.0000 

1.0000 

Table  4.2. 1-5.  Compression  Test  1-Reduced  Newton  Results. 


Number  of  Steps 

I 

h 

b 

1 

91.0000 

91.0000 

1.0000 

Table  4.2.2- 1.  Tension  Test  Parameters. 


/„  =  0.9  psi,  e„  =  0.94 

a, I  =  Cij;  =  033  =  -6.92308x  10  *  psi,  0,2  =  0,3  =  Ojj  =  0 
Ae,,  =  Ae,,  =  Ae,,  =  9.15692xl0■^  Ae,,  =  As,,  =  AS23  =  0 


Table  4.2.2-2.  Tension  “Exact”  Results. 


1 

h 

b 

-2.07692x10-^ 

0.5 

1 

Table  4.2.2-3.  Tension  Test  Trapezoidal  Results  before  /„  Fix- 


Number  of  Steps 

/ 

b 

1 

2.0965 

0.0010 

~0 

Table  4.2.2-4.  Tension  Test  Trapezoidal  Results  with  /„  Fix. 


Number  of  Steps 

/ 

h 

b 

1 

2.0965 

0.6408 

~0 

2 

-3.0007 

0.6413 

0.0493 

4 

-1.6303 

0.6161 

0.0872 

100 

-0.1629 

0.5048 

1.0000 

Table  4.2.2-5.  Tension  Test  Closest  Point  Results. 


Number  of  Steps 

/ 

h 

b 

1 

-0.1145 

0.4961 

1.0000 

Table  4.2.2-6.  Tension  Test  Reduced  Newton  Results. 


Number  of  Steps 

/ 

h 

b 

1 

-0.1145 

0.4961 

1.0000 

Table  4.2.3- 1.  Shear  Test  l-Parameters. 


0 


Figure  4.2.3- 1.  Trapezoidal  Shear  Test  1-Results. 


O 


Figure  4.2.3-2.  Closest  Point  Shear  Test  1-Results. 


O 


Figure  4.2. 3-3.  Reduced  Newton  Shear  Test  1-Results. 


Table  4.2.4- 1.  Shear  Test  2-Parameters. 


O 


Figure  4.2.4- 1.  Trapezoidal  Shear  Test  2-Results. 


O 


Figure  4.2.4-2.  Closest  Point  Shear  Test  2-Results. 


O 


Figure  4.2.4-3.  Reduced  Newton  Shear  Test  2-Results. 


Table  4.3. 1-1.  Compression  Test  2-Parameters. 


J 

/„  =  81  (psi),  =  0.94 

Ou  =  Ojj  =  033  =  -15  psi,  0,2  =  0,3  =  O23  =  0  _ 

A£,,'=  A£22  =  Ae33  =  -1.0xl0■^  A£,2  =  Ae,3  =  AEj,  =  0 


Table  4.3. 1-2.  Compression  Test  2-Trapezoidal  Results. 


Number  of  Steps 

I 

Jo 

b 

1 

48.8177 

82.5852 

2.3143 

2 

48.8186 

82.5874 

2.3143 

4 

48.8188 

82.5878 

2.3143 

100 

48.8188 

82.5879 

2.3143 

Table  4.3. 1-3.  Compression  Test  2-Closest  Point  Results. 


Number  of  Steps 

/ 

Jo 

b 

1 

48.6779 

82.7207 

2.3342 

2 

48.7452 

82.6572 

2.3245 

4 

48.7855 

82.6193 

2.3193 

100 

48.8173 

82.5894 

2.3145 

Table  4.3. 1-4.  Compression  Test  2-Reduced  Newton  Results. 


Number  of  Steps 

/ 

Jo 

b 

1 

48.7494 

82.6532 

2.3240 

,  2 

48.7877 

82.6172 

2.3186 

4 

48.8017 

82.6040 

2.3167 

100 

48.8182 

82.5885 

2.3144 

Table  4.3.2- 1.  Shear  Test  3-Parameters. 


/„  =  81  psi,  e„  =  0.94 

a,,  =  a,,  =  Ojj  =  -15  psi,  a,,  =  a, 3  -  0^3  =  0 
Ae|,=Ae„=Ae33=  0,  Ae,2  =  -1.0x10%  Ae,3=Ae23=  0 


Figure  4.3.2- 1.  Trapezoidal  Shear  Test  3-Results. 


Figure  4.3. 2-2.  Closest  Point  Shear  Test  3-Results. 


Figure  4. 3.2-2.  Reduced  Newton  Shear  Test  3-Results. 


Table  4.3.3- 1 .  Shear  Test  4-Parameters. 


/„  =  81  psi,  e„  =  0.94 

<^11  =  ^22  =  <733  =  -15  psi,  0,2  =  -2.0  psi,  0,3  =  023  =  0 
Aei,=Ae22=Ae33=  0,  Ae.j  =  -1.0xl0■^  Ae,3=Ae23=  0 


Figure  4.3.3- 1.  Trapezoidal  Shear  Test  4-Results. 


Figure  4.3. 3-2.  Closest  Point  Shear  Test  4-Results. 


Figure  4.3. 3-3.  Reduced  Newton  Shear  Test  4-Results. 


Table  4.4. 1-1.  Nearly  Elastic  Compression  and  Shear  Test-Parameters. 


/„  =  81  psi,  =  0.94 

=  CT,;  =  C33  =  -10  psi,  a,,  =  a,3  =  Oi,  =  0 

Ae,,=A£„=Ae33=  -1.0x10%  Ae,,  =  1.0x10 ^  Ae,3=Ae,3=  0 


Figure  4.4. 1-2.  Nearly  Elastic  Test-Relative  Times  and  Substeps. 


Table  4.4.3- 1.  Softening  Shear  Test-Parameters. 


Figure  4.4. 3-2.  Softening  Shear-Test  Relative  Times  and  Substeps. 


5.  Constitutive  Model  Implementation 

Constitutive  models  are  generally  designed  to  be  implemented  within  finite  element 
programs.  The  manner  in  which  these  programs  solve  the  global  nonlinear  problem 
dictates  the  requirements  of  the  constitutive  model.  For  this  study,  the  implementation  will 
be  in  a  program  which  uses  the  Newton-Raphson  method  for  nonlinear  solutions.  The 
DYSAC2  finite  element  computer  program  is  a  dynamic  analysis  code  for  soils  and  was 
written  at  the  University  of  California  at  Davis  [Muraleetharan,  et.  al.,  1991]. 

This  section  will  be  divided  into  two  parts: 

1)  discussion  of  the  issues  of  interfacing  a  constitutive  model  into  a  typical  finite 
element  program  that  uses  the  Newton-Raphson  method,  and 

2)  results  of  implementing  the  Bounding  Surface  model  for  clays  into  the  DYSAC2 
program. 

5,1  Standard  Interface 

The  development  of  a  “standard”  interface  for  constitutive  models  creates  a  conflict 
between  modular  programming  and  the  means  for  convenient  debugging  (for  an  analyst  to 
debug  a  “bad”  solution  as  opposed  to  the  debugging  required  for  implementing  a 
constitutive  model).  These  two  opposing  philosophies  dictate  the  level  of  information  that 
should  be  provided  through  the  constitutive  model  interface.  Modular  programming 
requires  the  minimum  amount  of  information  in  the  interface,  whereas  debugging  often 
requires  more. 

Modular  programming  is  based  on  a  “need-to-know”  philosophy  and  provides  only 
the  information  required  to  perform  its  function.  The  benefit  of  this  minimal  interface  is 
that  it  prevents  additional  information  (especially  in  weakly-typed  languages,  such  as 
FORTRAN)  from  being  corrupted.  It  also  provides  the  easiest  means  for  adding  new 
constitutive  models  to  existing  programs  because  the  user  is  not  required  to  track  the 
meaning  of  extraneous  variables. 


Convenient  debugging,  on  the  other  hand,  requires  a  great  deal  of  information  that 
is  never  used  within  the  constitutive  model.  The  information  is  needed  for  tracking  the 
cause  of  the  problem  in  the  event  of  a  numerical  failure  within  the  model.  When  a 
constitutive  model  fails  (e.g.,  does  not  converge)  the  user  can  extract  information  (such  as 
the  integration  point,  element,  iteration  and  time  step  numbers)  to  determine  whether  an 
input  parameter  was  specified  wrongly  or  simply  that  the  mesh  was  not  detailed  enough  in 
a  particular  area.  This  information,  however,  is  not  required  for  the  execution  of  the 
constitutive  model. 

In  order  to  provide  the  best  of  both  philosophies,  two  interfaces  are  proposed: 

1)  the  constitutive  branching  routine,  and 

2)  the  standard  constitutive  model. 

The  constitutive  branching  routine  is  generally  called  from  the  integration  point 
routine  and  contains  the  logic  for  deciding  which  material  set  belongs  to  the  element.  The 
constitutive  model  can  be  called  from  five  separate  areas  in  the  program  during  the  analysis 
phase  (depending  on  the  program’s  architecture): 

1)  reading  of  constitutive  properties, 

2)  initialization  of  internal  variables, 

3)  stresses  and  tangent  moduli  calculation  during  the  Newton-Raphson  iteration, 

4)  updating  internal  variables,  and 

5)  stress  reporting. 

Often  the  updating  is  merely  the  swapping  of  a  temporary  and  permanent  internal 
variable  array,  and  stress  reporting  is  accessing  the  stored  stresses.  At  this  level,  sufficient 
debugging  information  should  exist  to  track  a  numerical  problem  in  the  constitutive  array. 
The  branching  routine  is  proposed  to  have  both  debugging  and  constitutive  information 


with  an  interface  defined  as; 


branching  routine  (analysis  phase,  time  step  number,  iteration 
number,  element  number,  integration  point 
number,  constitutive  information) 

The  constitutive  properties  and  the  initial  values  of  the  internal  variables  are 
obtained  during  the  reading  phase.  This  step  is  done  once  for  each  separate  material  type 
and  often  doesn’t  go  through  the  branching  routine  (since  element  information  is  not 
required  for  debugging).  A  typical  read  routine  interface  is  defined  as: 

read  routine  (properties,  internal  variables) 

output: 

properties  =  constitutive  property  array 
internal  variables  =  internal  variable  array. 

The  constitutive  information  required  will  be  defined  by  the  initialization  and 
Newton-Raphson  iteration  requirements.  A  typical  initialization  routine  interface  is: 

initialization  routine  (properties,  initial  stresses,  internal  variables, 

external  variables) 


input: 

properties  =  constitutive  property  array 
initial  stresses  =  initial  stress  state  array 
external  variables  =  external  variable  array 


output: 


internal  variables  =  internal  variable  array. 

During  the  New^ton-Raphson  iteration,  the  constitutive  routine  should  produce  the 
change  in  stress,  tangent  moduli  and  the  change  in  internal  variables  given,  the  constitutive 
model  properties,  initial  stresses,  existing  strains,  internal  and  external  variables  and  the 
increment  in  strain.  The  routine  should  be  written  so  that  when  numerical  errors  occur  (or 
are  about  to  occur)  they  are  trapped  and  it  exits  to  the  branching  routine  where  the  error  and 
necessary  debugging  information  can  be  reported.  This  is  known  as  defensive 
programming  and  prevents  the  code  from  crashing  within  the  constitutive  routine  where  no 
debugging  information  is  available.  A  typical  constitutive  routine  interface  is  defined  as: 

constitutive  routine  (properties,  strains,  strain  increment, 

stresses,  internal  variables,  external  variables, 
tangent  moduli,  error,  convergence  information) 


input: 

properties  =  constitutive  property  array 
strains  =  existing  strain  array 
strain  increment  =  new  strain  increment  array 
input  and  output: 

stresses  =  previous  stress  state  array  (overwritten  with  the  new  stresses) 
internal  variables  =  internal  variable  array  (overwritten  with  new  array) 
external  variables  =  external  variable  array  (possibly  overwritten) 
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output: 


tangent  moduli  =  the  stress  -  strain  gradient, 


I  matrix 

n+1 


error  =  convergence  errors 

convergence  information  =  whether  or  not  convergence  occurred. 


The  external  variables  (e.g.,  temperature,  T)  are  often  input  only.  If  the  analysis  is 
coupled  with  these  variables,  the  main  program  might  also  require  additional  gradients 


associated  with  them  to  be  defined 


d{flux) 

de 


As  shown  above,  the  stress  and  internal  variable  arrays  are  overwritten.  Often  a 
temporary  array  is  maintained  for  both  the  stresses  and  internal  variables.  Once  the  global 
solution  has  converged,  then  the  updating  simply  consists  of  copying  the  temporary  array 
into  the  permanent  array.  If  the  temporary  arrays  need  to  be  maintained  at  the  global  level, 
separate  arrays  of  the  previous  and  new  values  can  be  specified  at  the  interface.  This 
allows  for  faster  global  updating  (i.e.,  a  single  loop  over  all  of  the  update  arrays  without 
going  down  to  the  element  level)  and  easier  restart  capabilities.  The  constitutive  routine 
interface  can  now  be  defined  as: 


constitutive  routine  (properties,  strains,  strain  increment, 

previous  stresses,  previous  internal  variables, 
previous  external  variables,  new  stresses, 
new  internal  variables,  new  external  variables 
tangent  moduli,  error,  convergence  information) 


133 


input: 

previous  stresses  =  previous  stress  state  array 

previous  internal  variables  =  previous  internal  variable  array 

previous  external  variables  =  previous  external  variable  array 

output: 

new  stresses  =  new  stress  state  array 

new  internal  variables  =  new  internal  variable  array 

new  external  variables  =  new  external  variable  array 

When  temporary  arrays  are  used,  stress  reporting  at  the  integration  points  becomes 
a  matter  of  simply  printing  the  stress  array.  Calculation  of  stresses  at  other  points  within 
the  element  (e.g.,  node  points)  is  generally  done  at  a  global/element  level  (i.e.,  curve 
fitting)  in  order  of  average  element  contributions  to  common  points.  A  graphical 
representation  of  this  interface  is  given  in  Figure  5.1-1.  The  dotted  lines  indicate  calls  that 
the  main  program  can  make  directly  without  going  through  the  branching  routine. 

5.2  Application  of  the  DYSAC2  Program  to  an  Embankment  Subjected  to 
an  Earthquake 

In  order  to  evaluate  the  performance  of  the  Reduced  Newton  version  of  the 
Bounding  Surface  clay  model  (Section  3),  the  model  was  implemented  into  the  DYSAC2 
finite  element  computer  program  [Muraleetharan,  et.  al.,  1991].  DYSAC2  is  a  dynamic 
soil  analysis  code  for  two-dimensional  plane  strain  problems  that  fully  couples  the 
governing  equations  of  a  saturated  porous  media  (two  phase).  The  implementation  was 
made  in  accordance  with  the  “standard  interface”  discussed  in  the  previous  section. 

The  problem  chosen  for  the  evaluation  of  the  numerical  model  involves  an  earthen 
clay  embankment  subjected  to  an  earthquake.  The  embankment  was  previously  analyzed 


using  DYSAC2  (with  base  shaking)  using  the  three-surface  clay  model  with  trapezoidal 
integration  and  comparing  the  results  to  the  centrifuge  results  [Muraleetharan,  et.  al., 

1994],  For  this  study,  the  single-surface  model  (ellipse)  was  used  with  the  properties 
given  in  Table  5.2-1  using  both  the  trapezoidal  and  Reduced  Newton  integrations.  A  more 
detailed  description  of  the  parameters  is  given  in  an  earlier  publication  [Herrmann  and 
Mish,  1983b].  The  interest  here  is  not  to  compare  the  results  to  the  centrifuge  test  results 
(although  a  direct  comparison  to  the  experimental  results  indicates  a  slight  improvement  in 
accuracy),  but  to  compare  the  performance  of  the  Reduced  Newton  integration  to  the 
trapezoidal  integration.  The  finite  element  dicretization  is  shown  in  Figure  5.2-1. 

The  centrifuge  model  was  spun  up  to  80g,  allowed  to  consolidate  and  then 
subjected  to  base  motion.  The  initial  stresses  (created  by  the  centrifuge  spin  up)  in  the 
finite  element  model  were  calculated  using  SAC2  [Herrmann  and  Mish,  1983c;  Herrmann 
and  Kaliakin.,  1987b].  Since  the  current  version  of  DYSAC2  does  not  have  this 
capability,  the  initial  stresses  were  graciously  provided  by  Dr.  Muraleetharan.  DYSAC2  is 
then  used  to  calculate  the  embankment  response  to  the  base  motion.  The  input  base 
acceleration  is  shown  in  Figure  5.2-2. 

The  comparison  of  the  embankment  predictions  using  the  trapezoidal  and  Reduced 
Newton  methods  can  be  evaluated  by  looking  at  various  nodal  and  element  responses.  The 
analysis  used  1,008  time  steps  with  a  time  increment  of  1.5625x10  '*  seconds.  An 
additional  analysis  was  performed  using  the  trapezoidal  method,  but  with  twice  the  number 
of  steps  (Trapezoidal  x2),  for  accuracy  comparisons.  The  horizontal  and  vertical 
accelerations  and  displacements  are  shown  in  Figures  5.2-3  and  5.2-4,  respectively.  The 
accelerations  of  all  three  analyses  match  extremely  well.  The  displacements  show  a  slight 
discrepancy  near  the  end  of  the  analysis,  with  the  Reduced  Newton  method  closely 
following  the  more  accurate  trapezoidal  run.  Excess  pore  water  pressure  histories  for 
selected  elements  are  shown  in  Figure  5.2-5,  and  show  a  similar  variation  with  the 
Reduced  Newton  being  more  accurate.  Contour  plots  of  the  excess  pore  water  pressure  for 


135 


both  methods  are  shown  in  Figure  5.2-6  and  indicate  little  difference  between  the  two 
methods.  The  times  as  recorded  on  a  PowerMac  8100/100  are  given  in  Table  5.2-2  and 
shows  the  Reduced  Newton  method  slightly  faster.  The  global  iterations  per  step  are 
shown  in  Figure  5.2-7  and  the  results  are  quite  similar. 

5.3  Application  of  the  DYSAC2  Program  to  an  Embankment  Subjected  to  a 
Shock 


Another  analysis  was  conducted  of  the  clay  embankment  with  a  shock  applied  at  the 
base.  The  same  material  properties  were  used  and  are  given  in  Table  5.2-1.  The  model 
was  analyzed  at  the  same  80g  level  (created  by  the  centrifuge)  so  that  the  same  initial 
conditions  could  be  used.  The  base  shock  is  shown  in  Figure  5.3-1 . 

The  analysis  used  5,000  time  steps  with  a  time  increment  of  5.0x10'^  seconds. 
Similar  to  the  earthquake  study,  an  additional  analysis  was  performed  using  the  trapezoidal 
method,  but  with  twice  the  number  of  steps  (Trapezoidal  x2),  for  accuracy  comparisons. 
The  horizontal  and  vertical  accelerations  and  displacements  are  shown  in  Figures  5.3-2  and 
5.3-3,  respectively.  The  accelerations  are  shown  for  the  first  0. 1  second  for  clarity  and  ail 
of  the  methods  match  extremely  well.  The  displacements  show  a  slight  divergence  near  the 
end  with  the  trapezoidal  and  Reduced  Newton  methods  drifting  slightly  from  the  more 
accurate  Trapezoidal  x2  analysis.  Excess  pore  water  pressure  histories  for  selected 
elements  are  shown  in  Figure  5.3-4  and  show  results  similar  to  the  displacements. 

Contour  plots  of  the  excess  pore  water  pressure  for  both  methods  are  shown  in  Figure  5.3- 
5  and  indicate  little  difference  between  the  methods.  The  times  as  recorded  on  a  PowerMac 
8100/100  are  given  in  Table  5.3-1  and  show  the  Reduced  Newton  method  to  be  slightly 


faster. 


Table  5.2-1.  Bounding  Surface  Model  Parameters. 


Parameter 

Description 

Value 

X 

Slope  of  virgin  compression  line 

0.25 

K 

Slope  of  swell/compression  line 

0.05 

K 

Slope  of  critical  state  (triaxial  space) 

0.88 

M^/M, 

Ratio  of  CSL  slopes  extension  to  compression  (3D  only) 

1.0 

V 

Poisson's  ratio  or  shear  modulus 

0.3 

li 

Transition  I  for  log  to  linear 

101.4  kPa 

P  atm 

Atmospheric  pressure 

30.4  kPa 

R 

Shape  of  bounding  surface 

2.4 

C 

Projection  center  parameter 

0.0 

Elastic  zone  parameter 

1.0 

m 

Exponent  for  shape  hardening  function 

0.02 

Hardening  parameter  associated  with  compression 

3.0 

H/H, 

Ratio  of  hardening  parameters,  ext/  comp  (3D  only) 

1.0 

Ho 

Hardening  parameter  associated  with  I  axis 

2.0 

a 

Parameter  controlling  magnitude  of  hardening 

1.2 

w 

Parameter  controlling  decrease  of  hardening 

5.0 

Table  5.2-2.  Embankment  Subjected  to  an  Earthquake  Timing  Results. 


Method 

Time  (sec) 

Time  (min) 

Trapezoidal 

3252 

54.2 

Reduced  Newton 

3223 

53.7 

Table  5.3-1.  Embankment  Subjected  to  a  Shock  Timing  Results. 


Method 

Time  (sec) 

Time  (min) 

Trapezoidal 

12052 

200.9 

Reduced  Newton 

11789 

196.5 

Figure  5.2-1.  Finite  Element  Discretization  of  Embankment. 
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Figure  5.2-2.  Input  Base  Acceleration  for  Earthquake. 
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Figure  5.2-3.  Comparison  of  Acceleration  Histories  of  Various  Nodes  Due  to  Earthquake. 
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Figure  5.2-5.  Comparison  of  Excess  Pore  Water  Pressure  Histories  of  Various  Nodes  Due 

to  Earthquake. 
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max  4.1338E+01 


Figure  5.2-6a.  Contours  of  Excess  Pore  Water  Pressure  Due  to  Earthquake,  Trapezoidal 

Method. 
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Figure  5.2-6b.  Contours  of  Excess  Pore  Water  Pressure  Due  to  Earthquake,  Reduced 

Newton  Method. 
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Figure  5.3-2.  Comparison  of  Acceleration  Histories  of  Various  Nodes  Due  to  Shock. 
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Figure  5.3-3.  Comparison  of  Displacement  Histories  of  Various  Nodes  Due  to  Shock. 
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Figure  5.3-4.  Comparison  of  Excess  Pore  Water  Pressure  Histories  of  Various  Nodes  Due 

to  Shock. 
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Figure  5.3-5a.  Contours  of  Excess  Pore  Water  Pressure  Due  to  Shock,  Trapezoidal 


Method. 
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Figure  5.3-5b.  Contours  of  Excess  Pore  Water  Pressure  Due  to  Shock,  Reduced  Newton 


Method. 


6.  Conclusions  and  Recommendations 


6.1  Conclusions 

The  Bounding  Surface  Plasticity  theory  provides  a  means  to  gradually  transition 
from  elastic  to  plastic  behavior  within  a  plasticity  setting.  Other  approaches  using  multi¬ 
surface  plasticity  have  been  proposed  to  achieve  the  same  behavior  [Mroz,  1967;  Krieg, 
1975;  Eisenberg  and  Phillips,  1971],  This  capability  is  especially  useful  in  cyclic  loading 
cases  where  the  change  in  volume  occurs  gradually  over  a  number  of  cycles.  Classical 
plasticity  models  consist  of  a  yield  surface  defined  in  stress  space  where  stress  paths  within 
the  surface  are  elastic  and  paths  on  the  surface  are  plastic.  Cyclic  behavior  that  occurs 
within  the  surface  is  elastic  and  can  never  generate  the  observed  plastic  behavior.  Unlike 
classical  plasticity  models,  the  surface  defined  in  a  Bounding  Surface  model  is  a  bound 
rather  than  a  yield  surface  and  therefore  allows  plasticity  to  occur  within  the  surface,  not 
just  when  the  stress  path  reaches  it.  This  feature  presents  difficulties  when  trying  to 
implement  the  model  using  conventional  numerical  techniques  which  were  developed  for 
yield  surface  models. 

This  study  investigates  numerical  techniques  for  specifically  implementing  the 
Bounding  Surface  model  for  clays.  The  first  goal  was  to  lay  out  the  Bounding  Surface 
theory  in  a  framework  that  would  allow  implementation  using  conventional  numerical 
methods.  Section  1.7  describes  an  approach  that  explicitly  defines  the  loading  surface  in 
terms  of  an  additional  internal  variable  so  that  the  model  fits  the  conventional  plasticity 
framework. 

The  Closest  Point  Projection  method  (Section  2)  is  one  of  the  modem  techniques 
for  numerically  evaluating  conventional  plasticity  models.  The  use  of  an  elastic  predictor 
and  a  plastic  corrector  makes  it  very  well  suited  for  yield  surface  plasticity  that  has  distinct 
elastic  and  plastic  zones.  Adding  an  additional  internal  variable  and  the  loading  surface  rate 
equation  allows  a  direct  application  of  this  method  to  the  Bounding  Surface  model.  Since 


the  method  explicitly  satisfies  the  consistency  condition,  it  often  provides  a  more  accurate 
answer  than  the  trapezoidal  method  (used  previously  with  the  Bounding  Surface  model  for 
clays)  in  fewer  steps.  However,  it  requires  solving  for  and  inverting  a  nonsymmetric 
matrix  (8x8  Jacobian)  in  every  iteration  at  the  material  point  level.  This  Jacobian  is 
eventually  used  to  constmct  the  consistent  tangent  moduli,  so  the  tangent  moduli  are 
effectively  being  calculated  in  every  iteration.  Also,  since  the  hardening  modulus  {K^) 
grows  rapidly  near  the  elastic  surface  or  projection  point,  the  loading  surface  rate  equation 
becomes  ill-behaved  and  requires  substepping  for  nearly  elastic  steps. 

The  Reduced  Newton  method  was  developed  from  some  of  the  concepts  in  the 
Closest  Point  Projection  method.  It  maps  the  stresses,  strains  and  consistency  condition 
into  invariant  space  and  employs  the  Newton-Raphson  method  to  solve  the  reduced  set  of 
equations.  The  loading  surface  relationship  developed  in  Section  1 .7  is  still  used,  but  in  a 
more  general  form  that  avoids  having  the  numerical  stability  problems.  The  method  also 
has  the  advantage  of  solving  for  the  stress  increment  first,  then  calculating  the  tangent 
moduli  only  for  the  converged  stress  solution. 

Both  the  Reduced  Newton  and  the  Closest  Point  Projection  methods  show  the  same 
behavior  in  terms  of  accuracy  (Section  4).  In  addition,  both  methods  tend  to  be  more 
accurate  than  the  trapezoidal  method  for  integration  of  large  plastic  strains,  because  they 
explicitly  include  the  consistency  condition.  For  robustness,  both  methods  (like  the 
trapezoidal  method)  employ  uniform  substepping  within  the  global  strain  increment,  which 
assures  a  consistent  level  of  accuracy  within  the  global  finite  element  mesh. 

The  Reduced  Newton  method  calculates  a  nearly  exact  Jacobian  (i.e.,  “consistent 
tangent  stiffness  matrix”),  while  the  trapezoidal  method  uses  a  tangent  stiffness 
approximation.  In  light  of  this  fact,  is  somewhat  surprising  that  the  two  methods  require 
substantially  the  same  number  of  global  iterations.  This  is  most  likely  due  to  the  fact  that 
the  dynamic  analysis  required  very  small  time  steps. 


The  development  of  the  Reduced  Newton  method  requires  extensive  mathematical 
manipulations  and  computer  coding.  The  method  is  best  suited  for  well  established  models 
for  which  reducing  computational  time  is  desired.  On  the  other  hand,  the  Closest  Point 
Projection  method  is  easier  to  implement  and  is  ideal  for  testing  constitutive  models  under 
development. 

From  both  an  accuracy  and  cost  standpoint,  the  trapezoidal  method  compared  more 
favorably  to  the  other  two  methods  than  expected.  As  a  consequence,  it  may  be  desirable 
to  attempt  to  correct  the  problem  of  straying  outside  the  bound  by  strictly  enforcing  the 
consistency  condition  within  the  trapezoidal  method. 

6.2  Recommendations  for  Future  Research 

This  study  has  suggested  a  number  of  possibilities  for  future  research  and 
improvements  on  the  Bounding  Surface  model  for  clays  and  its  numerical  implementation. 

Recommendations  for  additional  research  in  the  numerical  implementation  of 
Bounding  Surface  models  are  as  follows: 

a)  Investigate  a  trapezoidal  algorithm  that  explicitly  satisfies  the  consistency  condition. 

b)  Investigate  a  Closest  Point  Projection  algorithm  where  the  interval  is  divided  into  N 
substeps,  the  first  N-l  of  equal  length  and  the  last  with  half  of  the  spacing.  For  the 
first  N-\  substeps,  the  finite  difference  equations  would  be  written  at  the  center  of 
the  substep  (in  an  attempt  to  improve  accuracy),  while  for  the  last  substep  a 
backward  formula  would  be  used  so  consistency  at  the  end  point  is  exactly 
satisfied.  It  is  expected  that  this  method  might  demonstrate  the  favorable 
convergence  characteristics  of  the  trapezoidal  method. 

c)  Extend  the  Reduced  Newton  method  to  three  invariants. 

d)  Further  improve  robustness  of  the  methods  by  computationally  efficient  techniques, 
such  as,  variable  substepping. 


Specific  recommendations  for  additional  research  on  the  existing  clay  model  include 
the  following: 

e)  Determine  a  more  computationally  efficient  expression  for  the  hardening  modulus 
inside  the  bound.  The  modulus  described  in  Section  1.6  is  extremely  complicated 
and  requires  finite  difference  approximations  for  its  derivatives.  A  simpler, 
computationally  efficient  modulus  would  improve  the  speed  significantly.  This 
simpler  form  would  come  from  one  or  more  of  the  following  recommendations. 

f)  Define  a  bounding  shape  that  more  closely  follows  the  critical  state  line.  Part  of  the 
reason  for  the  complicated  hardening  modulus  is  that  the  elliptical  shape  extends  far 
beyond  the  critical  state  line.  Other  shapes  in  the  various  regions  have  been  tried 
with  some  success,  but  they  have  significantly  increased  both  the  complexity  and 
the  computational  time  of  the  model. 

g)  Recast  the  consolidation  law  to  eliminate  the  two  solution  regions  defined  by  the 
transitional  volumetric  stress  (I,). 

h)  Consider  a  model  without  a  purely  elastic  region.  This  type  of  model  would  be 
more  appropriate  for  clays  and  would  eliminate  the  two  separate  solution  regions. . 

i)  Consider  treating  the  unloading  implicitly  (as  was  done  for  the  Closest  Point 
algorithm)  instead  of  explicitly  (i.e.,  calculation  of  the  intersections).  It  is  assumed 
that  the  global  step  sizes  would  be  sufficiently  small  that  this  approximation  would 
be  valid  (which  is  often  the  case  in  earthquake  analyses). 

j)  Incorporate  inherent  and  induced  anisotropy  into  the  model. 
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Appendix  A.  Derivatives  for  Single  Ellipse  Bounding  Surface 


The  Bounding  Surface  function  for  a  single  ellipse  model  is  described  in  Section 
1.4  [Kaliakin,  1985;  Herrmann,  et.  al.,  1985;  Herrmann,  et.  al.,  1987].  The  function  is 
given  as: 


F  =  l^+{R-lf 
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As  described  in  Section  1.8,  the  Bounding  Surface  function  is  differentiated  with 
respect  to  the  “image”  stress  using  the  chain  rule  and  the  invariants.  This  can  be  written  as: 
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The  first  derivatives  are  given  as: 
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The  second  partial  derivative  of  F  with  respect  to  the  “image”  and  current  stresses  is 


given  as: 

d-F  d^F  dl  dF  d^l  d^F  dJ  dF  d^J 

-  —  - — - j-  - — - f - ;;;;; - - -  -f*  - - - 

<9rt.<9(7y  dld(7i^i  do^j  dl  do^jdo^i  dJda^,  d(7y  dJ  do-.da^ 

d^F  dfl  dF  d^H 

+ - £-4- - —JZ — 

d^d(J,i  dOy  d^  dGyda,, 

The  terms  for  the  second  derivative  are  found  to  be: 


d^F 

dido,, 


2b6, 


dU 


do,jda 


kl 


=  0 


d^F  dJ  2J  dN  dfl  ^ 

dJdo„  ~  1,  do,,  N  dll  do,,  j 

d^J  _  1  ds,^ _ s,j  dJ 

do,jdo„  2J  do,,  2J^  do,^ 


d-F 

diido,. 


2bJ  dJ 

^  do,, 


'  dN  dll 
^N)  dll  do,. 


d'H  _  dK  y  X  ^ 
do,jdo„  do,,  do. 


do^ 


kl 


(A-8) 

(A-9) 

(A-10) 

(A-11) 

(A-12) 

(A-13) 

(A-14) 


where 


161 


It  is  especially  important  that  the  above  definition  of  the  invariants  that  involve  /t  be 
used  to  form  the  second  derivatives.  In  the  previous  form  (see  Section  1.8),  the  term 
cos  (a)  appears  several  times  in  the  denominator  within  these  second  derivatives.  This 


creates  numeric  problems  as  cos  (oc)  approaches  zero. 
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Appendix  B.  Plastic  Algorithms  B  and  C 

Another  approach  for  the  plastic  algorithm  was  developed  originally  for  the  Cam- 
Clay  model  [Herrmann,  1997].  This  method  solved  a  single  nonlinear  differential  equation 
to  evaluate  the  material  behavior.  Because  of  numerical  problems  near  the  volumetric  stress 
axis  (i.e.,  7  =  0)  and  the  apex  of  the  ellipse  (/  =  maximum),  different  forms  of  the  equation 
were  used  in  these  two  regions.  The  first  section  of  this  appendix  discusses  the  two 
separate  equations  that  need  to  be  solved  when  using  this  approach  for  the  stress  point 
algorithm.  The  second  section  describes  the  calculation  of  the  consistent  tangent  moduli. 


Appendix  B.l  Stress  Point  Algorithm 


Section  3.4  develops  the  nonlinear  differential  equations  that  describe  the  material 
behavior.  The  volumetric  equations  are  given  as: 


-p7  )  =  A6I-4  — 


(/>/;)  (B.l-la) 


1  7, 


n+\ 


P  h 


(7<7,)  (B.l-lb) 


The  deviatoric  equation  is  given  as: 


+2G„„Af;^2  +  = 


l  +  2G„+]Ar— 


(B.1-2) 


The  above  equations  are  sufficient  for  the  Cam-Clay  model.  For  the  Bounding 
Surface  plasticity  model,  an  additional  equation  is  required  to  describe  the  plasticity 
behavior  within  the  surface.  The  equation  used  is: 
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(B.1-3) 


The  current  algorithm  described  in  Section  3.4  solves  for  I  and  J  using  the  two 
nonlinear  Equations  B.1-1  and  B.1-2  and  then  Equation  B.  1-3  to  solve  for  the  plasticity 

parameter  {'f}.  The  original  implementation,  which  is  consistent  with  the  Cam-Clay  model 

implementation  [Herrmann,  1997],  solved  Equation  B.1-3  and  either  Equation  B.1-1  or 
B.  1-2  for  /  and  b.  The  equation  not  used  in  this  process  was  used  to  find  the  plasticity 
parameter.  The  choice  of  Equation  B.  1-1  or  B.  1-2  is  based  upon  numerical  considerations 
and  is  described  later  in  this  section. 

In  both  cases,  the  first  equation  is  Equation  B.1-3  and  gives  the  Newton-Raphson 
residual: 
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Algorithm  B  uses  Equation  B.1-2  for  the  second  equation.  Equation  B.1-1  is  then 
used  to  solve  for  the  plasticity  parameter.  Residual  2  is  defined  as: 
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A0--^ 


P  1 1 
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Numerical  problems  arise  in  this  algorithm  near  the  apex  of  the  Bounding  Surface 
(J  =  maximum)  because  the  denominator  of  Equation  B.1-6  approaches  zero.  To  resolve 
this  problem.  Algorithm  C  exchanges  the  second  residual  and  the  plasticity  parameter 


equations.  Residual  2  and  the  plasticity  parameter  for  Algorithm  C  are  defined  as: 


(/>/;)  ,  (B.l-7a) 
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(/</,)  (B.l-7b) 


(B.1-8) 


The  plus  sign  before  the  radical  is  required  to  assure  a  positive  plasticity  parameter. 
Numerical  problems  arise  in  this  algorithm,  however,  near  the  volumetric  axis  (i.e.,  J  =  0). 
Solving  for  the  plasticity  parameter  using  Equation  B.1-8  results  in  inaccurate  estimates 
because  of  the  division  by  small  values  of  J  (the  f  divisor  within  the  radical). 

Because  of  the  nature  of  these  equations,  the  solution  realm  is  broken  into  two 

zones: 

1)  Algorithm  B  is  used  for  solutions  in  the  bottom  half  of  the  Bounding  Surface 
(near  the  hydrostatic  axis),  and 

2)  Algorithm  C  is  used  for  solutions  in  the  top  half  of  the  Bounding  Surface  (near 
the  apex). 

The  dividing  point  for  the  algorithms  is  arbitrarily  taken  as  the  half  point  and  is 
based  upon  the  second  invariant  value  at  the  beginning  of  the  step: 


Algorithm B:  JI<^N^ 
2 


Algorithm  C:  -^b 


(B.1-9) 


M 

where  N  =  —j=  (angle  of  the  critical  state  line  in  I-J  space) 

V27 


M  =  angle  of  the  critical  state  line  in  p-q  space. 
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Similar  to  Section  3.4,  the  independent  variables  /  and  b  are  chosen.  The  Newton- 
Raphson  method  requires  derivatives  of  the  residuals  with  respect  to  the  independent 
variables.  The  derivatives  of  the  first  residual  (Equation  B.1-4)  with  respect  to  and 
are  given  as: 
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The  derivatives  of  the  second  residual  (Equation  B.1-5)  and  the  plasticity  parameter 
(Equation  B.  1-6)  for  Algorithm  B  are  given  as: 
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dy 
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The  derivatives  of  the  second  residual  (Equation  B.1-7)  and  the  plasticity  parameter 
(Equation  B.1-8)  for  Algorithm  C  are  given  as: 
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The  convergence  criterion  chosen  for  this  study  is  that  the  change  in  residuals  must 
be  smaller  than  a  user-defined  percentage  of  absolute  value  of  the  total  strain  increment. 
The  first  residual  is  in  terms  of  the  plasticity  parameter.  In  order  to  compare  its  value  at  the 
magnitude  of  strain,  the  residual  and  its  derivatives  are  multiphed  by  the  following: 
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Algorithm  C  uses  the  volumetric  Equation  B.  1-7  for  the  second  nonlinear  equation, 
which  is  already  in  terms  of  strain.  Algorithm  B,  however,  uses  the  deviatoric  Equation 
B.  1-5,  which  is  in  terms  of  stress.  Its  residual  is  scaled  to  the  magnitude  of  strain  by 
dividing  it  and  its  derivatives  by  the  initial  shear  modulus.  The  new  residual  and 
derivatives  are  given  as: 


diC,  _  1  ^  ^  1  dR^ 
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(Algorithm  B)  (B .  1  -2 1 ) 


The  norm  as  described  in  Section  3.4  is  based  upon  a  percentage  of  the  input  strains 


and  is  defined  as; 
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where  |/?||  =  norm  on  the  residuals. 


The  Newton-Raphson  algorithm  is  the  same  as  described  in  Section  3.4  and  is 
given  as: 
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where 
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<9/?2 

dR^ 

dl 

db 

Residual  2  in  the  Jacobian  ('R)  is  dependent  on  the  algorithm  being  used.  The  algorithm  is 


shown  in  Box  B.1-1. 


Box  B.1-1.  Algorithms  B  and  C  for  Plastic  Steps. 


1.  Set  Substep  Level:  N  =  1 

2.  Substep  Level  Loop:  m  = 

ASy 

Ae,..  = — - 

‘j.  ^ 

3.  Initialize  to  Beginning  of  Step: 

K+I  ~  h>  ^n+l  ~ 

4.  Substep  Loop:  k  =  1,N 

5.  Determine  Volumetric  Form: 

il>li)  Linear-Linear  Log-Linear 

6.  Determine  Algorithm: 

jI  <  —N—  Algorithm  B  N—  Algorithm  C 

2  b,  2  b, 

7.  Iteration  Loop:  it  =  i'it^ 

8.  Calculate  Residuals  and  Local  Jacobian:  R,, 

9.  Solve  for  Increment  in  Variables:  A/„^/,  Ab^,^, 

10.  Calculate  and  Test  Norm:  IF  (\\R\\<  Toler,) 

TRUE:  Go  to  step  1 1 

FALSE:  Increment  variables,  b^^,.  Go  to  step  7 

11.  Calculate  Stresses  and  Check  Substeps: 

JF  {k<N)  TRUE:  Go  to  step  4  FALSE:  Go  to  step  12 

12.  Compare  stresses  with  previous  substep  level: 

K,.. or  [">•"»«] 

TRUE:  EXIT  FALSE:  N  =  2N,  Go  to  step  2 


Appendix  B.2  Plastic  Contribution  to  the  Global  Jacobian 


Once  the  stress  algorithm  has  converged,  the  plastic  contribution 


da. 

- —  to  the 


global  Jacobian  matrix  must  be  calculated.  The  local  Jacobian  will  be  evaluated  in  terms  of 
the  deviatoric  and  volumetric  components,  similar  to  the  elastic  Jacobian.  The  volumetric 
component  can  be  obtained  by  noting  that  the  residuals  are  a  function  of  the  independent 
variables  (I,  b)  which,  in  turn,  are  a  function  of  the  strain  increment.  The  derivatives  can 
be  obtained  by  an  application  of  the  chain  rule  (i.e.,  holding  the  independent  variables 
constant  and  then  adding  the  derivatives  with  respect  to  the  independent  variables).  At 
convergence,  the  residual  is  equal  to  zero  and  the  derivatives  can  be  written  as: 
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The  local  Jacobian  ('P)  is  the  same  as  that  required  in  the  local  Newton-Raphson 
iteration  (B.  1-23).  The  derivatives  of  the  independent  variables  with  respect  to  the  strain 
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increments  can  be  found  by  solving  Equation  B.2-2: 
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The  derivatives  of  the  residuals  with  respect  to  the  strain  increments  are  calculated 
holding  the  independent  variables  at  the  end  of  the  step  b^^,)  constant.  It  is  important 
to  note  that  the  shear  modulus  is  a  function  only  of  the  current  volumetric  stress 

and  therefore  is  also  fixed.  For  both  algorithms,  the  derivative  of  the  first  residual  is 
given  as: 
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The  derivatives  of  the  second  residual  for  Algorithm  B  and  that  of  the  plasticity 
parameter  are  given  as: 
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For  Algorithm  C,  the  derivatives  of  the  second  residual  and  the  plasticity  parameter 


are  given  as: 
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The  derivatives  in  Equations  B.2-5  through  B.2-8  are  as  follows: 
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The  derivative  of  the  deviatoric  stress  with  respect  to  the  increment  in  strain 
in  Equation  B.2-16  is  evaluated  at  the  beginning  of  the  step  or  substep.  If  the 


initial  portion  of  the  step  was  elastic,  then  this  derivative  is  calculated  using  Equation  3.3-4. 


If,  on  the  other  hand,  the  initial  portion  of  the  step  was  plastic,  then  the  derivative  comes 
from  the  previous  plastic  step.  The  deviatoric  stress  derivative  at  the  end  of  the  step  is 
defined  as: 
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This  derivative  requires  the  derivatives 


31 


and 


,  which  are  given  by 


Equation  B.2-3. 


Appendix  C.  Strain  Calculations  for  “Exact”  Solutions 


Appendix  C.  1  Volumetric  Problems 

In  order  to  verify  that  the  numerical  implementation  of  the  Bounding  Surface 
models  are  calculating  the  correct  results,  several  “exact”  solutions  were  constmcted.  Exact 
is  in  quotes  because  in  some  cases  (such  as  the  volumetric  tests)  the  exact  solution  can  be 
found.  In  general,  however,  numerical  integration  is  required.  The  “exact”  solutions  are 
constructed  by  starting  with  a  stress  point  on  the  Bounding  Surface,  giving  the  surface  and 
the  stress  point  a  known  movement  together  and  calculating  the  strains  required  to  achieve 
the  movement.  Since  the  stress  point  starts  on  the  surface,  this  solution  could  also  be  used 
to  test  a  Cam-Clay  model. 

The  “exact”  solution  is  constructed  by  assuming  a  linear  expansion  with  time  (?)  of 
the  Bounding  Surface  size  {Iq)'. 

=  (C.1-1) 


where  /„  =  initial  Bounding  Surface  size. 

The  derivative  of  this  Bounding  Surface  expression  with  respect  to  time  is: 


/  =1  (C.1-2) 

() 

The  previously  defined  rate  equation  for  the  Bounding  Surface  hardening  [Kaliakin, 
1985]  is  given  as: 


(C.1-3) 


where  Vo  =  I +eo  =  initial  specific  volume 


Co  =  initial  void  ratio 
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X  =  normal  consolidation  line 


K=  elastic  load-reload  line 


//  =  Limit  stress  for  linear  relationship 


7=  plasticity  parameter  (loading  index). 


Consider  the  case  when  the  value  of  ^  is  less  than  the  transition  stress  (//).  The 
Bounding  Surface  relationship  reduces  to: 


i  y  ^ 


X-K 


dl 


(C.1-4) 


where  I  =  b[l  - 1^)  + 


b  =  measure  of  distance  between  stress  point  and  surface 


L  =  Cl 


C  =  material  constant  defining  the  projection  center  location. 
The  single  ellipse  Bounding  Surface  function  is  defined  [Kaliakin,  1985]  as: 


F  =  I^+{R  +  lf 


ay 

UJ 


2  ,  -  2-R  ,2 

R^"  ^  R  " 


(C.1-5) 


The  intercepts  of  the  Bounding  Surface  with  the  volumetric  axis  (7  axis)  are: 


(C.l-6) 


The  derivative  of  the  Bounding  Surface  function  (Equation  1.4-1)  with  I  is: 


^  =  21-^1. 
dl  R  " 


(C.1-7) 
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Substituting  Equation  C.1-7  into  Equation  C.1-4  and  recalling  Equation  C.1-2: 


1(0  = 


3  h 

X-K 


r 


f27-|/„(0 

V 


\ 


=  1 


(C.l-8) 


Equation  C.l-8  can  now  be  solved  for  the  plastic  parameter  {f): 


7  = 


V3  hj 


1 


n-\m 


(C.1-9) 


The  rate  equation  for  the  plastic  volumetric  strain  (6^)  is  given  by  the  volumetric 
form  of  Equation  1.2-2  [Kaliakin,  1985]: 


r  =  3rf 


(C.1-10) 


Substituting  in  the  expressions  for  the  plastic  parameter  and  derivative  (Equations 
C.1-7  and  C.1-9)  yields: 


d”  = 


X-K 


(C.1-11) 


The  plastic  volumetric  strain  can  be  obtained  by  integrating  the  rate  with  respect  to  f. 


=  +  c 


(C.1-12) 


where  C  =  constant  of  integration. 


The  constant  of  integration  is  evaluated  by  noting  that  ^  =  0  at  t  =  0,  and  therefore 
C  =  0.  Thus  the  equation  for  the  plastic  volumetric  strain  is: 


d’’{t)  = 


X-K^ 
V  h  J 


(C.1-13) 


179 


The  total  volumetric  strain  consists  of  both  the  plastic  and  elastic  strain.  This 
relationship  is  given  by  the  familiar  decomposition  of  strains  written  in  volumetric  form: 

d{t)  =  G\t)  +  e'^{t)  (C.1-14) 

In  order  to  determine  the  elastic  volumetric  strain,  the  stress  history  is  required. 

For  volumetric  stresses  (/)  less  than  the  transition  stress  (//),  the  bulk  modulus  {K)  is 
constant: 

K  =  ^ 

(C.1-15) 

The  volumetric  stress  rate  is  defined  as: 

i{i)  =  3Ke‘(t)  (c,i.i6) 

Equation  C.1-16  can  be  integrated  for  the  elastic  volumetric  strain  (ff).  This  is  a 

straightforward  integration  since  the  bulk  modulus  is  constant.  The  constant  of  integration 

is  zero  since  /  =  0  at  0^  =  0.  Substituting  the  elastic  strain,  the  plastic  volumetric  strain  (0^ 

in  Equation  C.1-13)  and  the  bulk  modulus  (K  in  Equation  C.1-15)  into  Equation  C.1-14, 
the  total  strain  can  be  written  as: 

e(,)  =  EJKjtzJi]  , 

K'^oh)  (C.1-17) 

Consider  a  volumetric  tension  test.  The  stress  point  originates  on  the  volumetric 
axis  intercept  of  the  Bounding  Surface  on  the  tension  side  (for  /?  >  0)  and  remains  there 
throughout  the  test.  I(t)  can  be  defined  by  the  tension  side  intercept  of  the  bounding 
surface  (first  term  of  Equation  C.1-6).  Substituting  this  into  Equation  C.1-17  results  in  the 


volumetric  strain  as  a  function  of  the  Bounding  Surface  size  only: 


X-K 


'/  J 


(C.1-18) 


This  equation  for  strain  highlights  a  problem  with  the  Bounding  Surface  rate 
relationship  given  in  Equation  C.1-3.  When  considering  the  case  where  the  Bounding 
Surface  collapses  to  zero  (i.e.,  /„  ->  0  ),  t  cannot  be  less  than  -  /„  .  Therefore  the  total 

tensile  volumetric  strain  in  this  limit  cannot  exceed: 


/  \ 

1 

II 

1 

V  /  1 

1  ^0  h  J 

(C.1-19) 


Tensile  strains  exceeding  this  value  will  cause  the  model  to  have  numeric 
difficulties.  In  order  to  avoid  this  problem,  one  approach  is  to  reconsider  the  Bounding 
Surface  rate  (Equation  C.1-3).  Since  there  was  no  experimental  evidence  cited  for  the 
transition  from  a  log  to  linear  rate  form,  the  linear  portion  is  dropped.  The  rate  equation  for 
lo  (previously  defined  in  Equation  C.1-3)  is  now  defined: 


3  v„ 


X-K 


7 


dF 

dl 


(C.1-20) 


Note  that  this  form  of  the  rate  equation  is  valid  on  either  side  of  the  transition  stress. 
The  plastic  parameter  is  now  solved  for: 


7  = 


^X-K^ 


3v. 


1 


(C.l-21) 


Substituting  the  plasticity  parameter  into  the  rate  equation  for  the  plastic  volumetric 
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strain  (6^)  (Equation  C.1-10)  results  in: 


0”  = 


V  J 


1 


4(0 


(C.1-22) 


Integrating  to  get  the  plastic  volumetric  strain  gives: 


0^{t)  = 


\  ) 


In  (4(0)  +  C 


(C.1-23) 


The  constant  of  integration  is  obtained,  noting  that  ^  =  0  at  r  =  0: 


C  =  - 


V  J 


In  (4(0) 


(C.1-24) 


The  equation  for  the  plastic  volumetric  strain  as  a  function  of  time  is: 

e\t)  = 


r 

In 

C  r  / 

4(0 

J 

V 

(C.l-25) 


Substimting  this  relationship  into  the  equation  for  the  total  volumetric  strain 
(Equation  C.  1-14)  for  tensile  tests  results  in: 


0{t)  = 


K 


in 


\  '^<’1  J 


(C.1-26) 


When  considering  the  case  where  the  Bounding  Surface  collapses  to  zero 
(i.e.,  /„  — >  0),  t  =  -  /„  .  The  total  tensile  volumetric  strain  is  now  given  as; 

e(-4,)  = 


(C.1-27) 
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Now  in  order  to  collapse  the  Bounding  Surface  to  zero,  an  infinite  value  of  tensile 
volumetric  strain  is  required.  Although  the  precision  of  the  machine  will  limit  this  strain, 
the  formulation  will  avoid  the  numerical  difficulties  of  the  previous  formula. 

For  volumetric  stresses  (7)  greater  than  the  transition  stress  (7/ ),  the  bulk 
modulus(70  is  now  a  function  of  the  volumetric  stress: 

3/s:  (C.1-28) 

The  integration  of  the  elastic  volumetric  strain  (Equation  C.1-16)  results  in  a 
slightly  more  complex  expression: 


(C.1-29) 


(C.1-30) 


Consider  a  volumetric  compression  test.  The  volumetric  stress  (7)  would  originate 
on  the  compression  side,  intercept  with  the  Bounding  Surface  and  remain  on  the  surface  as 
the  bound  expands.  Substituting  for  the  volumetric  stress  (i.e.,  second  term  of  Equation 
C.1-6),  the  total  strain  becomes  a  function  of  the  Bounding  Surface  size: 


(C.1-31) 

The  volumetric  compression  test  in  Section  4.2.1  used  Equation  C.1-30  to  calculate 
the  strains.  The  Bounding  Surface  parameters  are  given  in  Table  C.1-1.  The  void  ratio 
was  given  as  =  0.94.  The  problem  definitions  are  given  in  Table  C.  1-2. 

The  volumetric  tension  test  in  Section  4.2.2  used  Equation  C.1-26  to  calculate  the 
strains.  The  Bounding  Surface  parameters  and  the  void  ratio  are  the  same  as  the 


6{t)  =  —  In 
t/„ 


I  •¥  t 


compression  test.  The  parameters  defining  the  problem  definitions  for  the  tension  tests  are 


given  in  Table  C.1-3. 


Table  C.1-1.  Bounding  Surface  Model  Parameters 


Traditional  Model 

Surface  Configuration 

Parameters 

Parameters 

II 

P 

R  =  2.6 

K  =  0.05 

M,=  1.05 

V  =  0.2 

Table  C.1-2.  Volumetric  Compression  Problem  Definitions 


Bound  Size, 

Volumetric  Stress,  I 

Volumetric  Strain,  d 

Beginning 

81 

81 

0 

Final 

91 

91 

-0.00840075 

Table  C.1-3.  Volumetric  Tension  Problem  Definitions 


Bound  Size,  /„ 

Volumetric  Stress,  I 

Volumetric  Strain,  6 

Beginning 

0.9 

0.207692 

0 

Final 

0.5 

0.115385 

0.0274708 
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Appendix  C.2  Shear  Problem 

A  problem  is  designed  to  test  the  behavior  of  the  model  in  shear.  In  order  to  make 
finding  the  exact  solutions  tractable,  some  simplifying  assumptions  are  made: 

1)  the  first  invariant  remains  constant  (i.e.,  =  /„),  and 

2)  the  third  invariant  is  zero  (5  =  0). 

Similar  to  the  volumetric  cases,  a  linear  variation  in  lo  is  assumed: 


(C.2-1) 


As  derived  in  the  previous  section  (Equation  C.1-21),  the  plastic  parameter  {f)  is 


given  as: 


y  = 


J 


1 


K 


(C.2-2) 


The  plastic  volumetric  strain  (Equation  C.1-25)  is  given: 


e^it)  = 


V  y 


In 


V  J) 


(C.2-3) 


Since  the  volumetric  stress  (/)  is  assumed  to  remain  constant,  the  total  volumetric 
strain  can  be  expressed: 


m  =  d^it) 


(C.2-4) 


The  second  image  stress  invariant  ( 7)  is  obtained  as  a  function  of  time  by  recalling 
the  equation  for  the  Bounding  Surface  (Equation  C.1-5)  and  noting  that  F  =  0: 


7(0  = 


^  A-P+-  /„(r)/-— ^ 

(F-l)V  R 


nit) 


(C.2-5) 
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The  equation  for  N  is  defined  as: 


N  =  - 


2N 


N, 


1- 


c  V 


N. 


3S(S 


u 


(C.2-6) 


Noting  that  the  third  invariant  is  assumed  to  be  zero  (i.e.,  S  =  0),  the  equation  for  N 
becomes  constant  and  can  be  written  as: 


N  = 


2 


!  +  • 


N. 


(C.2-7) 


For  tests  originating  and  remaining  on  the  surface,  the  first  and  second  invariants 
can  be  simplified  to: 


I  =  I  J  =  J 


(C.2-8) 


In  order  to  ensure  that  1  =  constant  and  5  =  0,  the  following  assumption  on  the 
deviatoric  stresses  is  made: 


■^11  ~  ^22  ~  '^33  “  '^13  ~  ^23  ~  ^  ^ 


The  shear  stress  ((Tj2)  can  be  written  as: 


(C.2-9) 


^12(0  —  ‘^12(0  ~  *^(0 


The  rate  equations  for  the  plastic  deviatoric  strains  are  defined  as: 


e^  =  r^ 


(C.2-10) 


(C.2-11) 
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The  derivatives  of  the  Bounding  Surface  function  (F  in  Equation  C.1-5)  with 
respect  to  the  deviatoric  stresses  are  given  as; 

dF  _dF  dJ  ^  dF  d^i 

dJ  ds,j  (C.2-12) 

where  ^  “  sin(3cK) 

Specific  derivatives  of  the  Bounding  Surface  function  are  given  as: 


2j_ 


(C.2-13) 


djJ. 


=  -2 


2^^ 

dfx 


(C.2-14) 


dN 

dF- 


_ _ V 


N. 


c  / 


IN^ 


(C.2-15) 


dF 


1-- 


V 


N. 


N 


(C.2-16) 


The  derivatives  of  the  invariants  with  respect  to  the  deviatoric  stresses  are: 


dJ  dJ  dJ  dJ  dJ 


ds^  j  ds'^  ds-^'^  ds'^ 


=  0 


'23 


(C.2-17) 


dJ  \_ 
2 


ds 


12 


(C.2-18) 


dF  _  dF 

ds22  27 


(C.2-19) 
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dfl 

ds- 


33 


A 

j 


(C.2-20) 


dfl  ^  dfl  ^  d/i  ^  ^ 
dS^'y  ^*^23 


(C.2-21) 


Upon  combining  these  terms,  the  derivatives  of  the  Bounding  Surface  function  with 
respect  to  the  deviatoric  stresses  are  defined  as: 


^5*1  j  (55*22 


2  _  g 

V 

N  K 


J{t) 


(C.2-22) 


dF 


ds- 


=  V3(i?-l) 


33 


\-!F' 

f  V  N,j 
N 


J{t) 


(C.2-23) 


dF  1 


<95, 


=  2(ff-i)"^y(<) 


12 


(C.2-24) 


dF  dF 


^5,3  ds 


=  0 


23 


(C.2-25) 


Substituting  in  the  plastic  parameter  (Equation  C.2-2),  the  second  invariant 
(Equation  C.2-5),  and  the  derivatives  (Equations  C.2-22  through  C.2-25)  into  the 
deviatoric  strain  rate  (Equation  C.2-1 1),  the  rates  can  be  written  as: 


pP  -  gP  = 
^22 


12 


1-K 

V  y 


K 


1- 


2  ,  r  r  A  r  ^  ^  72 


N.. 


m 


(C.2-26) 


e’’  =-2  e’’ 

^33  z,  Cij 


(C.2-21) 
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A.-K 


V  y 


N 


(C.2-28) 


To  obtain  the  plastic  deviatoric  strains,  these  rates  must  be  integrated  thus  requiring 
numerical  evaluation  of  the  last  term.  With  this  approximation,  the  deviatoric  plastic  strains 
can  be  expressed  as  a  function  of  time.  The  total  deviatoric  strains  are  obtained  by  recalling 
that: 


(C.2-29) 

The  shear  modulus  (G)  is  a  function  of  the  bulk  modulus  {K)  which,  in  turn,  is  a 
function  of  the  first  stress  invariant  (/).  Assuming  that  I  =  constant,  the  shear  modulus  is 
also  constant.  From  Equation  C.2-9,  the  total  deviatoric  strains  are  obtained: 

=  ef,  ^22  =  4  ^33  =  <3  (C.2-30) 


=  ^  +  (C.2-31) 

Finally  the  total  strain  histories  are  calculated  using  Equation  C.2-4: 


(no  sum) 


(C.2-32) 


=  (C.2-33) 

The  Shear  Test  1  in  Section  4.2.3  used  the  Bounding  Surface  parameters  given  in 
Table  C.1-1 .  The  void  ratio  was  given  as  =  0.94  and  the  volumetric  stress  was  set  to  /  = 
81  psi.  The  parameters  defining  the  problem  are  given  in  Table  C.2-1. 

Shear  Test  2  in  Section  4.2.4  also  uses  the  Bounding  Surface  parameters  given  in 
Table  C.1-1.  The  void  ratio  was  given  as  =  0.94,  the  volumetric  stress  was  set  to 


/=  60  psi  and  the  initial  shear  stress  was  given  as  (Jj2  =  3.63731  psi.  The  parameter 
definitions  are  given  in  Table  C.2-2. 


Table  C.2-1.  Shear  Test  1  Definitions 


Bound  Size,  /„ 

Volumetric  Strain,  6 

Shear  Strain,  e,2 

Beginning 

81 

0 

0 

Final 

91 

0.00540048 

-0.0114129 

Table  C.2-2.  Shear  Test  2  Definitions 


Bound  Size,  /„ 

Volumetric  Strain,  6 

Shear  Strain,  £,2 

Beginning 

81 

0 

0 

Final 

91 

0.00715132 

-0.0177091 

